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SUMMARY 


The  study  is  a  part  of  an  effort  being  directed  toward  solving  the 
problem  of  a  composite  material  subjected  to  general  oblique  loading. 

This  analysis  was  conducted  in  order  to  find  the  internal  micro¬ 
mechanics  of  a  fiber-reinforced  composite  due  to  transverse  normal  loading. 
Special  emphasis  has  been  given  to  studying  the  stress  distribution  near 
free  surfaces,  which  led  to  solving  a  three-dimensional  elasticity  problem. 
The  numerical  method  of  finite  elements  has  been  employed  in  this  analysis. 
On  the  other  hand,  it  was  necessary  to  study  the  behavior  far  from  free 
surfaces.  For  this  purpose,  a  two-dimensional  program  was  used.  Findings 
from  these  two  approaches  were  anticipated,  showing  that  stress  conditions 
become  two  dimensional  a  relatively  short  distance  from  the  end  of  the 
composite.  Extensive  parametric  studies  have  been  performed  from  the 
combined  outputs  of  these  two  schemes.  Significant  diagrams  exhibiting 
elastic  properties  of  different  composite  materials  have  been  obtained. 
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INTRODUCTION 


Only  in  the  ideal  case  is  the  loading  of  a  reinforced  composite  mate¬ 
rial  in  the  directions  of  the  reinforcements.  Normally,  oblique  loading 
is  encountered  by  the  composite.  To  obtain  the  effect  of  oblique  loading, 
we  must  combine  axial,  shear,  and  transverse  loadings.  By  doing  so,  we 
can  establish  the  stresses  in  the  reinforcements  and  in  the  matrix  as  they 
actually  occur;  we  can  then  use  these  stresses  for  the  development  of  a 
strength  theory  based  on  such  reliable  data. 

This  report  presents  the  case  of  transverse  loading  of  a  uni¬ 
directional  composite.  The  two  different  regions  considered  are  (1)  the 
region  close  to  the  end  of  the  fibers,  where  three-dimensional  treatment 
is  applied,  and  (2)  regions  far  from  the  end  of  the  fibers,  where  the  end 
perturbations  are  not  more  effective  and  where  a  plane  solution  is 
applied . 

In  both  cases,  the  finite  element  method  was  applied  successfully. 

For  the  plane  problem,  we  also  intended  to  find  a  solution  by  using  the 
so-called  point-matching  or  collocation  method  which  was  applied  in 
the  solution  of  the  longitudinal  shear  problem.  However,  the  convergence 
of  this  method  was  very  poor  for  this  case  of  transverse  loading. 

Appendix  I  gives  the  equations  utilized,  along  with  a  description  of  the 
form  in  which  the  point-matching  method  was  used. 

The  plane  problem  was  solved  in  both  plane  stress  and  plane  strain 
conditions  by  using  the  finite  elemert  method.  Displacements,  stresses, 
and  transverse  moduli  for  different  types  of  composites  are  given  in 
this  report. 


THREE-DIMENSIONAL  SOLUTION 


The  numerical  approach  used  to  evaluate  a  specified  three-dimensional 
boundary  problem  of  linear  elasticity  represents  a  first-order  approxima¬ 
tion  method.  This  method  is  known  as  the  method  of  finite  elements.  Its 
concept  is  based  on  the  assumption  that  the  stress  within  small-volume 
elements  of  the  body  is  constant.  This  is  precisely  true  if  we  consider 
infinitesimal  volume  elements.  For  small-volume  elements,  this  assumption 
of  constant  stress  proves  to  be  a  good  model  representation  of  reality. 

It  is  convenient  to  consider  small  tetrahedrons  as  volume  elements 
because  three-dimensional  space  can  be  subdivided  into  sets  of  tetra¬ 
hedrons  in  a  simple  way.  The  assumption  of  constant  stress  within  each 
elementary  tetrahedron  is  equivalent  to  the  assumption  of  linear 
displacement-vector-distibution  within  the  elementary  tetrahedron. 

Compatibility  and  equilibrium  conditions  introduced  at  the  tetra¬ 
hedron  node  points  lead  to  a  system  of  linear  algebraic  equations  whose 
solutions  represent  displacements  at  the  node  points  of  the  tetrahedrons. 
From  these  displacements,  we  can  determine  the  stresses  within  each 
tetrahedron. 

THE  STIFFNESS  MATRIX  OF  A  TETRAHEDRON  ELEMENT 


In  order  to  obtain  the  proper  working  equations,  we  have  to  consider 
an  elementary  tetrahedron  first  and  find  the  pertinent  relations  as  far 
as  stress,  nodal  forces,  ard  nodal  displacements  are  concerned. 

Figure  1  shows  a  general  elementary  tetrahedron  with  the  nodes  N;  , 
Na  ,  N3  ,  and  N*  .  We  assume  that  the  three  coordinates  §i  >  'Hi  »  and 
Qi  (i  “  1,2, 3, A)  are  known  for  each  node  point  N^  (i  =  1,2, 3, 4)  . 


Figure  1.  The  Elementary  Tetrahedron. 
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Again,  we  can  represent:  the  displacement  vector  by  a  linear  vector 
function  of  the  local  coordinates  ?  ,  T)  ,  and  £  .  The  origin  of  the 
local  coordinate  system  coincides  with  node  Ni  .  In  this  fashion,  we 
put  the  three  displacement  components  at  point  P(§,T).C)  >  as  follows: 


u  =  a.  + 

x  1 

aa  ?  + 

a3 

H  + 

a4  C 

(1) 

u  =  a5  + 

y  5 

ae  l  + 

a7 

H  + 

aS  C 

(2) 

c 

N 

II 

+ 

aio!  + 

al  1 

H  + 

ai  sC 

(3) 

The  coefficients  ai>****ai2  are 

constants 

and 

subject 

to  variance  with 

the  geometric  configuration  of  the  tetrahedron  as  well  as  the  displace¬ 
ment  configuration  at  the  nodes.  Let  ux.  ,  Uy.  ,  and  uz^  (i  *■  1,2, 3,4) 
represent  displacement  components  of  the  four  node  points.  Then, we  obtain 
the  following  12  relations  from  equations  (1),  (2),  and  (3): 


u  «=  a:  +  aa  ^  +  a3  Hi  +  a4  £i 

xi 

U  =  al  +  as  ^2  +  a3  ^2  +  a4  Ql 

*2 

llx,  =  f1  +  a2  ^3  +  a3  ^3  +  a4  Cs 

=  ^  +  3j  ?4  +  a4  "f  C4 

uyi  =  ae  +  ae  ?i  +  37  Hi  +  ae  Cl 

u  =  a,.  +  ai  Hj  +  as  Cj 

y8  & 

“y3  ■  S  +  a.  53  +  a,  \  *  a,  C, 

•  S  +  H  5<  +  ^  \  +  %  c. 

U7  =  a0  +  aio?!  +  ^i1!-  +  ^aCi 

zi 

uza  =  a9  +  a10?a  +  +  ^  3^8 

uZ:j  =  S  +  +  ajjCa 

U  =  a9  +  310^4  +  a! !  ^4  +  auk  (4) 

Z4 


3 


In  matrix  form,  relations  (4)  can  be  expressed  as  follows: 


1  00000000000 

1  F2  Tls  Ca  0  0  0  0  0  0  0  0 

1  ?a  \  C3  00000000 

1  U  ^4  C»  0  0  0  0  0  0  0  0 

OOOOluOOOOOO 
0  0  0  0  ]  ?,  T)a  Ca  0  0  0  0 

0  0  0  0  1  §3  Tj3  q3  0  0  0  0 

OOOOl^T^^OOOG 
000000001  000 
0  0  0  0  0  0  0  0  1  5a  T|s  £a 

0  0  0  0  0  0  0  0  1  §3  £3 

0  0000000  1  5*  k 


Equation  (5)  can  be  written  symbolically,  such  as 

ju|  =  [A]  | a |  (b) 

where  {u}  are  the  12  displacements  at  the  four  nodes,  {a}  the  {l2  x  l} 
matrix  of  the  configuration  coefficients  a  ,...,a  ,  and  [A]  the  { 1 2  x  12} 
matrix  as  demonstrated  in  equation  (5).  1 

The  six  strain  components  ex  »  ey  »  ®z  •  €xy  >  exz  >  anc^  eyz  arc 
obtained  immediately  from  equations  (1),  (2),  and  (3),  as  follows: 


“  ^7 
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_*  +  — * 
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xz 
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du  9u 
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(7) 


We  denote  {e}  as  the  (6  X  l)  matrix  with  the  elements  ex  ,  ev  ,  e 


-xy 


€xz  ,  and  Cyz  and  write  equation  (7)  symbolically,  as  follows: 
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(8) 


If  we  use  the  short  form  symbol  (Dj  for  the  (6  X  12)  matrix  in  equa¬ 
tion  (8),  we  can  write  equation  (8)  in  the  symbolic  short  form: 


M  -  MM 


(9) 


Finally,  we  can  obtain  the  six  stress  tensor  components  ax  ,  Cy  ,  az  , 
aXy  »  cxz  ,  ayZ  from  Hooke's  law: 
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with  the  meanings  of  X  and  G  , 


and 


2  (1+v) 

(ID 

(1+v)  (1  -2v) 

(12) 

where  E  is  the  modulus  of  elasticity  and  v  is  Poisson's  ratio. 

The  symbol  [C]  stands  for  the  {6  X  6}  matrix  of  equation  (10), 
which  is  the  elasticity  matrix.  Equation  (10)  therefore  is,  in  short 
symbolic  form, 

M  ■  [C]  |  =  i  (O) 

Substitution  of  {e}  by  means  of  equation  (9)  follows: 

ja|  -  [C][D]  |  a  |  (14) 

From  equation  (6),  we  can  express  {a}  in  terms  of 

jaj  -  [A]'1  juj  (15) 


where  [a]"1  represents  the  inverse  matrix  of  [a]  .  Therefore,  equation 
(14)  is 


jaj  -  [cMa]-1  juj 


(16) 


Expression  (16)  is  the  working  equation  for  calculating  the  stress 
tensor  {a}  .  In  order  to  do  this,  it  is  necessary  to  Inow  all  12  dis¬ 
placement  components  at  the  4  nodes.  In  general,  whether  or  not  each  one 
of  them  is  given  input  quantities  depends  on  the  kind  of  boundary  condi¬ 
tions  present.  As  a  matter  of  fact,  we  are  in  general  not  free  to  choose 
a  displacement  component  arbitrarily  when  the  corresponding  force  compo¬ 
nent  has  been  fixed.  Castigliano' s  theorem  expresses  this  proposition 
clearly,  stating  that  in  any  linear  elastic  system,  the  force  component 
acting  on  the  system  in  a  certain  direction  is  equal  to  the  first  deriva¬ 
tive  of  the  strain  energy  function  with  respect  to  the  displacement  in 
the  same  direction.  If  we  define  W  as  the  strain  energy  function  of 
the  system,  we  can  obtain  a  set  of  equations  which  relates  the  node 
force  component  array  to  the  displacement  components,  as  follows: 


6 


(17) 


W  proves  to  be  a  quadratic  form  in  fu},  and  therefore  it  is  the  right-hand 
side  of  a  linear  function  in  {u}  representing  the  desired  relationship 
between  forces  and  displacements.  Next,  we  have  to  obtain  the  strain 
energy  function  W  . 


The  strain  energy  of  an  elastic  body  is  the  volume  integral  over  the 
elastic  potential,  such  as 


2  j( 


a  e  + 
x  x 


o  e  + 
y  y 


a  e  + 
z  z 


a  e 
xy  xy 
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xz  xz 
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yz  yz 


dV  (18) 


T 

The  integrand  can  be  written  symbolically  [e]  {a}  ,  and  since  e  and  O 

are  constant  within  the  region  of  integration,  we  get 


where 


X3-X1 

ya-yi 

Z3-Zl 

xa-xx 

y3-yi 

Z3  -  z: 

X4  -xi 

y*  -yi 

Z4  _Z1 

£a  %  Ca 

§3  Ik  C3 

1) *  U 


(19) 


(20) 


is  the  volume  of  our  elementary  tetrahedron.  Keeping  in  mind  that  the 
volume  must  be  positive  regardless  of  the  selected  sequence  of  the  node 
points,  the  absolute  value  of  the  triple  product  in  equation  (20)  is 
indicated. 

Relation  (19)  is  the  desired  function,  and  all  that  is  left  is  the 
substitution  of  the  proper  linear  expressions  in  terms  of  {u}  ;  namely, 
expressions  (16)  and  (9)  in  conjunction  with  (15).  Substitution  of  {a] 
by  (19)  issues  the  modified  equation  (9). 

jcj  =  [D]  [A]*1  juj  (21) 

The  transposed  linear  matrix  {e}^  follows  after  applying  the  transposi¬ 
tion  rule  twice,  as  follows: 
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(22) 


HT  ■  {[*](*]■'■  M}1  -  {w-  H}T[»]T 

.  |u|T(tA]-‘)  *  Of 

With  equations  (22)  and  (16),  the  strain  energy  function  W  can  be 
written  in  the  following  quadratic  form: 

W  -  iv|u;T([A]-‘)T  [d]T  [C]  [D]  or  juj  (23) 


Equation  (23)  leads  to  the  desired  force  displacement  relationship  through 
the  second  Castigliano  theorem  (17).  Applying  the  differentiation  rules 
yields 

1*1  •  f(wl)  TmT [c][d][*]-‘  m 


+  !  HT(W‘)  ’M'MMW-1  »*> 


Repeated  application  of  the  transposition  rule  and  observation  that 
[C]^=  [C  ]  (the  elastic  matrix  is  a  symmetric  one)  shows  that  the  last 
right-hand  term  is  equal  to  the  first  term.  We, therefore,  get  the  system 
of  linear  equations 

|P|  -  v([A]f[D]T[C][D]  [A]'1  juj  (25) 

which  contains  the  set  of  working  equations  to  obtain  the  unknown  displace¬ 
ments  for  computing  the  stress  tensor  according  to  equation  (16)  .  We  call 

(K)  =  v([A]-»)  T  [D]T  [C]  [D]  [A]’1  (26) 

the  stiffness  matrix.  The  stiffness  matrix  [K]  is  a  symmetric  one  since 
the  matrix  of  elastic  constants  [C]  is  symnetric,  as  pointed  out  above. 

The  matrix  [A]  depends  only  on  the  local  coordinates  of  the  tetra¬ 
hedron  nodes.  Matrix  (26)  therefore  depends  only  on  the  local  coordinates 
of  the  tetrahedron  nodes  and  the  elastic  constants  of  the  elastic  medium 
within  the  tetrahedron. 

The  physical  significance  of  each  K-matrix  element  can  be  demonstrated 
by  considering  special  load  conditions  at  the  nodes  of  the  tetrahedron. 
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Let  us  assume  that  we  want  to  interpret  the  meaning  of  the  ij^ 
element  of  the  K-matrix  on  the  element  in  the  i*-*1  row  and  jt*1  column.  If 
we  set,  in  equation  (25),  all  displacement  components  but  the  one  equal 
to  zero  and  the.  displacement  component  equal  to  one  (single  displacement 
condition),  then  the  12  force  components  become  equal  to  the  12  coeff.  ients 
in  the  jt"  column.  The  value  of  the  coefficients  in  the  i*-*1  row  and  j1*1 
column  is  therefore  the  i*-^1  force  component  which  must  be  present  to  maintain 
all  displacement  components  at  zero  except  the  j*-*1  displacement  component, 
which  must  be  kept  at  unity.  In  general,  we  must  apply  12  force  components 
to  maintain  this  special  strain  condition  of  the  tetrahedron  in  correspond¬ 
ing  with  the  12  coefficients  in  the  column. 

The  j1-*1  force  in  particular  acts  in  the  same  direction  as  the  jth 
unity  displacement.  Sine  the  force  is  doing  work  on  the  sysi em,  the 
orientation  of  force  and  tisplacement  in  the  jth  direction  must  corre¬ 
spond.  This  means  that  all  diagonal  elements  of  the  K-matrix  must  be 
positive . 

In  many  cases,  it  is  the  single  unity  displacement  condition  that 
helps  us  to  visualize  intuitively  the  physical  significance  of  certain 
matrix  elements. 


THE  COMPOUND  STIFFNESS  MATRIX  OF  A  UTTICE 
FORMED  BY  A  SYSTEM  OF  TETRAHEDRONS 


Any  space  region  can  be  subdivided  into  systems  of  parallelepipeds 
and  each  parallelepiped  into  six  tetrahedrons.  In  this  way,  the  space 
region  can  be  built  up  by  a  system  of  tetrahedrons  in  the  same  way  that  a 
two-dimensional  region  can  be  covered  by  a  net  of  triangles. 

If  we  consider  an  interior  tetrahedron  —  that  is,  one  completely 
surrounded  by  four  other  tetrahedrons  and  whose  node  points  are  not 
boundary  points  of  the  space  regions  —  then  it  can  be  shown  that  each 
node  point  must  be  common  with  node  points  of  17  or  23  other  tetrahedrons 
in  the  neighborhood  of  the  tetrahedron  being  considered.  This  means  that 
the  forces  generated  at  a  lattice  point  of  interior  space  are  the  sum  of 
all  forces  generated  by  the  18  or  24  tetrahedrons, each  having  one  of  its 
node  points  coinciding  with  the  lattice  point  being  considered. 

Consequently,  we  can  now  see  that,  in  general,  the  force  at  a  lattice 
point  will  be  affected  by  the  displacements  at  each  node  of  each  of  the 
18  or  24  connected  tetrahedrons.  This  is  evident  since  all  proper  force 
equations  from  each  of  the  18  or  24  contributing  tetrahedrons  must  be 
summed  together  in  order  to  get  the  resultant  force  at  the  lattice  point. 
In  this  it  uner,  we  get,  in  general,  a  linear  relationship  between  the 
compounds  force  af  che  lattice  point  and  all  pertinent  displacement 
component.  _ 'pearing  in  all  the  18  or  24  tetrahedrons  connected  at  the 
lattice  point. 
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Let  us  suppose  that  che  i^1  lattice  point  is  connected  to  p  tetra¬ 
hedron,  having  the  numerical  sequence  tx  (i) , t£ (i), , . .  , tp (i )  .  The  elements 
of  this  sequence  are  functions  of  the  lattice  point,  as  indicated  by  the 
argument.  The  force  vector  at  the  ith  lattice  point  is  represented  by  its 
three  components 


p<‘> 

,w 

rs 

p(i) 

r3 


(27) 


The  stiffness  matrix  of  the  tetrahedron  t  shall  be  expressed  by 


(28) 


and  the  line  elements  of  this  matrix  belonging  to  the  fqrce  components 
1,2,3  at  the  lattice  point  i  are 


kM,I 

r  ( t ,  i )] 

k<c'IJ] 


(29) 


The  compound 


force  at  the 


ith 


lattice  point  is  therefore 


where 
the  tj 


is  the  column  nodal 
tetrahedron. 


displacement  vector  corresponding  to 


The 'total  matrix  that  is  obtained  in  this  fashion  for  all  N  lattice 
points  is  denoted  by  [K* ]  ,  which  is 


p 

M  -  £ 

[^(tjU),!)' 

j“l 

[^(tjd),^ 

(31) 


where  i  =  1,2,3,. ,N  . 
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It  is  obvious  that  the  configuration  of  tetrahedrons  connected  to  a 
lattice  point  i  will  be  different  when  i  is  a  boundary  point  itself. 

The  column  vector  jtj^j  will  contain  lesser  elements  than  the  corres¬ 
ponding  matrix  for  interior  lattice  points. 


We  have  distinguished  between  the  force  components  belonging  to  a 
lattice  point  i  and  in  a  certain  direction  (1,  2,  or  3).  This  scheme  has 
the  advantage  of  allowing  us  to  know  with  which  element  we  are  dealing.  In 
order  to  use  one  index  number  to  identify  a  force  (or  displacement)  component 
but  still  contain  the  information  of  the  lattice  number  i  and  the  type  of 
direction  (J  =  1,  2,  3),  we  introduce  the  index  number  scheme: 

n  «  3(i  -  1)  +  j 


(i  =  1.2...N  ;  j  =  1,2,3)  .  In  this  way,  the  indices  of  components  form 
a  sequence  of  natural  numbers,  and  each  number  still  contains  the  informa¬ 
tion  of  being  an  x,  y,  or  z  component  as  well  as  the  lattice  index 
n  mber  i  .  This  information  can  be  extracted  easily,  since 


j  =  mod(n,3)  +  1 

where  mod(n,3)  is  the  remainder  of  division,  n  divided  by  3,  and 

i  =  +  1  if  mod(n,3)  ^  0  / 

oj.  ? 

i  =  if  mod(n,3)  =  0  J 

where  means  the  largest  integer  contained  in  the  quotient  ^  . 


With  matrix  (31),  the  system  of  force  components  acting  upon  the 
given  configuration  of  lattice  points  can  be  expressed  symbolically  by 


M  =  i>]  H 


Since  we  use  indexing  of  forces  and  displacements,  the  system  of  equations 
(33)  can  be  written  in  analytical  form,  such  as 


>*  =  Y.  K*.  u* 

1  fi  ^  J 
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where  i  “  l,2,3...M  .  We  assume  now  that  there  is  M  , 
of:  force  components,  as  well  as  displacement  components. 

will  be  one  of  M  x  M  . 


the  total  number 
The  matrix  [K*l 


Because  of  applied  forces  and  displacements  introduced  to  the  system, 
there  are  some  force  components  known  and  the  remaining  force  components 
unknown.  The  same  is  true  for  the  displacement  components  in  a  complemen¬ 
tary  sense.  If  there  are  M,  unknown  force  components,  then  there  are 
M  -  1^  known  force  components  and  there  will  be  M  -  unknown  displace¬ 

ment  components  while  there  are  K  known  displacements.  Let  us  arrange 
the  equations  (14)  such  that  the  first  equations  express  Lhe  equilib¬ 

rium  of  the  unknown  forces  and  the  rest  express  the  equilibrium  of  the 
known  forces.  In  the  same  fashion,  we  arrange  the  terms  in  each  equation 
such  that  the  first  terms  contain  known  d i splacements  and  the 

remaining  terms  contain  the  unknown  displacements.  Furthermore,  we 
distinguish  by  superindexing  the  known  and  unknown  displacements  and 
forces  with  (k)  and  (u)  .  Then  we  get  the  following  two  systems  of 
equations : 


E 


i-H+i 


(u) 
u  . 

J 


(35) 


where  i  *  1,2.. .M^  ,  and 


(k)  4 
u  .  + 

1 


n 

E 

=rr+i 


K*. 

ij 


(36) 


where  i»F\  +  l,fy+2,..M. 

The  last  equation  is  system  of  linear  equations  for  the  M  -  Mj 
unknown  displacements  u(u)  ,  where  i  -  H  +  1  ,  +  2 , . .M  . 


With  the  solution  of  equation  (36),  the  unknown  forces  can  be  computed 
from  equation  (35).  Equation  (36)  is  a  non homogeneous  linear  system  of 
(M  -  H)  equations: 


M 

£  **. 

(u) 

U  .  " 

J 

p(k>-  E  **.  '!k> 

1  U  > 

(37) 

where  i»M1  +  l,M1  +  2, 

.  .M  . 
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In  order  to  apply  numerical  methods  to  solve  this  system  of  equations, 
it  is  necessary  to  re-index  the  unknowns  by  consecutive  natural  numbers. 

This  can  be  easily  done  by  generating  a  number  sequence  (MSK  matrix)  that 
contains  the  indices  of  the  unknown  displacements  in  a  consecutive  order. 

The  arguments  of  this  sequence  are  now  consecutive,  natural  numbers  which 
are  the  working  indices  for  solving  equations  (37).  By  means  of  the  MSK 
matrix,  it  is  now  possible  to  identify  each  solution  element  with  the 
original  index.  The  main  purpose  of  this  scheme  is  to  provide  the  capa¬ 
bility  to  extract  the  information  of  location  (lattice  point  index)  and 
direction  (x,  y,  or  z)  of  the  displacement  from  the  equation  working  index. 
This  capability  is  required  in  many  instances;  e.g.,  to  satisfy  the  polar 
symmetry  conditions  and  to  select  the  proper  sequence  of  the  12  displacement 
elements  for  each  tetrahedron  in  order  to  compute  the  stresses.  Stress 
within  a  tetrahedron  is  computed  according  to  equation  (16).  The  stress 
distribution  for  the  three-dimensional  region  of  elastic  medium  is  then 
obtained  by  applying  certain  boundary  conditions. 

In  the  following,  the  boundary  problem  will  be  specified  and  the 
scheme  of  the  numerical  approach  described. 


DESCRIPTION  OF  THE  NUMERICAL  SCHEME 
Definition  .of  the  Boundary  Problem 

The  main  objective  of  this  study  is  to  investigate  the  effect  of  free 
surfaces  on  the  stress  distribution  within  a  fiber  matrix  composite  under 
loads  transverse  to  the  fiber. 


Figure  2  shows  this  con¬ 
figuration  schematically.  The 
stress  distribution  problem  here 
is  three  dimensional.  From  this 
figure,  we  can  observe  some 
symmetry  conditions  if  we  assume 
that  the  fibers  form  a  periodic 
network  of  hexagonals  with  one 
fiber  axis  intersecting  at  each 
apex  and  centroid  of  the  hexa¬ 
gonals  (hexagonal  A,B,C,D,E,F  is 
one  of  them) . 

Around  each  fiber  cross 
section,  there  is  a  hexagonal 
such  that  another  network  of 
hexagonals  is  formed.  This 
hexagonal  is  called  the  basic 
fiber-resin  element  of  the 
system.  In  Figure  2,  the 
hexagonal  A' , B' ,C' ,D' ,E '  ,F  ’ 
is  such  an  element. 


Figure  2.  Basic  Fiber-Resin 
Element. 
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Figure  3  depicts  the  top  view  of  the  intersection  of  the  fiber-resin 

element . 

We  consider  now  an  external  force  acting  in  a  plane  normal  to  line 
B',E'  but  at  a  distance  which  is  a  high  multiple  of  distance  B ' , E *  away 
from  M  .  We  call  this  the  transverse  load  upon  one  system.  Because  of 
symmetry  of  load  and  configuration  conditions,  the  displacements  at  points 
within  the  fiber  axial-parallel  planes  going  through  lines  Mj  ,  Mg  and 
,  P8  are  constant  and  in  the  direction  of  the  external  force.  They 
are  oriented  opposite  each  other  but  are  of  equal  magnitude.  All  points 
of  the  plane  normal  to  line  and  containing  point  Pj  must  have 

displacements  along  the  plane  they  generate.  The  same  displacement  con¬ 
ditions  exist  for  points  in  the  plane  parallel  to  the  one  described  above 
but  containing  point  M  .  Both  planes  are  planes  of  symmetry  for  the 
displacement  contribution.  In  this  way,  there  is  a  periodic  repetition  of 
the  displacement  pattern  which  exists  within  the  prism  above  the  rectangular 
cell  M,N,H  ,Pj  in  the  direction  of  the  four  planes  of  symmetry  which  are 
the  side  planes  of  the  prism.  Therefore,  we  can  study  this  transverse  load 
problem  by  the  following  boundary  conditions,  which  are  based  on  the  new 
coordinate  convention  for  the  basic  prism  shown  in  Figure  4. 


CTz(x,y,o) 

■ 

0 

(38) 

CTxz(x,y,o) 

a 

0 

(39) 

a  (x,y,o) 

- 

0 

(40) 

ux(±c,y,z) 

- 

±k 

(41) 

uy(x,±|,z) 

■ 

0 

(42) 

du 

-%r  <x»y»0 

- 

0 

(43) 

3u 

(x.y.t) 

— *■ 

0 

(44) 

uz(x,y  ,0 

— 

0 

(45) 

The  last  three  equations  indicate  that  the  condition  of  plane  strain 
will  be  approached  with  increasing  distance  from  the  free  surface  z  ■  0. 
The  distance  l  can  be  chosen  freely  such  that  it  will  be  more  than  four 
times  the  dimension  of  the  diagonal  of  the  base  triangle.  This  assumption 
has  been  based  mainly  on  St.  Venant's  principle.  Our  numerical  results 
verify  the  correct  selection  of  t  in  this  respect. 
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Within  the  region,  two  types  of  elastic  materials  are  considered. 

The  cylindrical  bodies  have  the  modulus  of  elasticity  Ef  and  Poisson's 
ratio  vf  .  Outside  of  the  cylindrical  domains,  the  associated  elastic 
constants  are  and  vm  . 

With  these  premises,  it  is  now  possible  to  compute  the  stresses  by 
applying  the  method  of  finite  elements.  In  order  to  do  this,  it  is  first 
necessary  to  subdivide  the  interior  of  the  above-specified  three-dimensional 
region  into  a  system  of  elementary  tetrahedrons. 

The  System  of  Elementary  Tetrahedrons 

In  Figure  4,  we  refer  to  the  free-boundary  plane  z  -  0  as  "ground 
floor."  Furthermore,  we  define  eight  equidistant  planes  parallel  to  the 
ground  floor,  and  we  call  the  nc^  plane  the  "nth  floor."  The  top  plane, 
z  «  (,  ,  is  referred  to  as  the  "tenth  floor."  These  10  planes  subdivide 
the  parallelepiped  into  nine  1  lyers,  each  having  the  height  dz  -  {,/ 9  . 

Next,  we  have  to  subdivide  each  layer  into  prisms  with  triangular 
bases,  since  each  prism  can  be  again  subdivided  into  three  tetrahedrons. 

The  simplest  way  to  do  this  is  to  subdivide  each  layer  in  the  same  way, 
which  means  that  we  erect  a  Bet  of  prisms  on  the  ground  floor,  each  prism 
having  the  height  l  .  Each  layer  will  thus  have  the  same  layout  of  tri¬ 
angular  bases  as  the  one  on  the  ground  floor.  Figure  5  shows  the  way  the 
ground  floor  has  been  subdivided  into  96  triangle  bases  (8  rows  with  12 
triangles  in  each  row).  Each  triangle  represents  a  prism  with  the  height 
of  the  layer  thickness  d2  .  The  top  triangle  of  the  prism  is  a  replica 
of  the  base  triangle  and  by  itself  is  the  base  triangle  for  the  next  prism 
above  the  second  layer.  In  this  fashion,  all  prisms  in  each  layer  above  a 
base  triangle  are  identical  geometrically.  The  base  triangles  are  numbered 
on  the  ground  floor  in  a  consecutive  manner.  Because  of  central  synmetry 
properties,  only  55  triangles  need  to  be  considered.  With  the  same  scheme, 
the  lattice  points  at  the  ground  floor  are  numbered  consecutively  in  each 
row  and  column.  Because  of  central  symmetry  properties,  the  last  point 
on  the  ground  floor  is  the  origin,  which  is  lattice  point  number  32. 
Continuation  of  the  numbering  starts  with  the  lattice  point  above  number  1 
in  the  same  manner  as  in  the  ground  plane  below.  In  this  way,  we  can 
calculate  the  corresponding  lattice  point  number  J  in  the  p*"  floor  above 
the  I1*1  lattice  point  by 


J  -  I  +  32  (p  -  1) 


(46) 


As  pointed  out  earlier,  the  indices  scheme  of  the  force  and  dis¬ 
placement  components  can  be  linked  with  the  indexing  scheme  of  the 
lattice  points.  Consequently,  the  qth  component  (q  -  1  is  x-component, 
q  ■  2  is  y-component,  and  q  ■  3  is  z-component)  of  the  vector  quantity 
has  considered  the  following  index  Q  : 

Q  -  3(J  -  1)  +  9  or  Q  -  l[l  +  32 (p  -  1)  -  l]  +  q  (47) 
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We  call  the  prism  above  the  base  triangle  between  two  consecutive 
floors  the  elementary  prism.  Since  a  triangular  prism  can  be  subdivided 
into  three  tetrahedrons  in  several  ways,  we  must  establish  some  basic 
rules.  Three  tetrahedrons  in  any  triangular  prism  are  generated  through 
intersecting  the  prism  by  two  planes.  These  two  planes  are  generated  by 
the  three  lines  of  diagonals  in  each  of  the  three  rectangles  forming  the 
mantle  of  the  prism.  See  Figure  6.  Furthermore,  it  is  necessary  that 
one  diagonal  line  intersect  the  other  two  diagonals.  In  this  way,  all 
three  diagonals  form  a  linked  train  of  lines  (train  of  diagonals) 
with  disconnected,  or  open,  ends.  The  first  basic  rule  we  apply  in  order 
to  generate  the  tetrahedrons  is  the  following:  The  diagonals  of  coincid¬ 
ing  rectangles  belonging  to  adjacent  prisms  coincide. 


Figure  6.  Tetrahedron  Forming  a 
Triangular  Prism. 
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Ac  this  point,  we  distinguish  each  point  of  the  triangle  base  as  to 
the  topological  property  of  the  train  of  diagonals  with  respect  to  the 
relative  position  to  the  triangle  base  proper.  A  lattice  point  of  a  base 
triangle  at  which  two  diagonal  lines  intersect  each  other  is  said  to  be 
of  a  topological  a-property.  When  two  diagonal  lines  intersect  at  the 
lattice  point,  then  the  lattice  point  considered  has  a  topological 
3-property.  Finally,  there  is  a  lattice  point  topological  y-property 
when  one  diagonal  goes  through  the  lattice  point  under  consideration 
and  the  other  diagonal  goes  through  the  lattice  point  above  or  below 
the  one  under  consideration. 

The  first  conclusion  from  this  is  that  a  base  triangle  must  contain 
points  of  all  three  topological  properties.  The  second  conclusion  that 
follows  from  this  topological  consideration  is  that  all  apexes  of  base 
triangles  having  one  common  lattice  point  are  of  the  same  topological 
type. 


According  to  the  last  theorem,  it  is  therefore  appropriate  to  think 
in  terms  of  topological  properties  of  lattice  points  only.  If  we  know 
the  topological  property  of  each  lattice  point  in  a  base  triangle,  we 
know  the  lattice  points  of  the  four  apexes  for  each  of  the  three  tetra¬ 
hedrons  above  the  base  triangle.  With  reference  to  Figure  6,  we  can 
state  the  following: 

Tetrahedron  I.  The  three  base  apexes  are  the  three  lattice 
points  of  the  base  itself.  The  fourth  apex  is  the  8-P°int 
in  the  next  floor  (3'  in  Figure  6). 

Tetrahedron  II.  The  three  base  apexes  are  the  a-point  at 
the  base,  the  S‘P°int  above  the  base  (8*  in  Figure  6),  and 
the  y-point  above  the  base  (y'  in  Figure  6).  The  fourth 
apex  is  the  y-point  in  the  base  triangle. 

Tetrahedron  III.  The  three  base  apexes  are  the  three  lattice 
points  of  the  triangle  above  the  base.  The  fourth  apex  is 
the  y-point  in  the  base  triangle. 

The  last  three  properties  contain  the  algorithm  to  generate  the  four 
lattice  numbers  corresponding  to  the  four  apexes  of  each  tetrahedron.  It 
is  therefore  necessary  to  know  the  topological  properties  of  each  lattice 
point.  In  the  numerical  program,  the  coordinates  of  each  ground  floor 
lattice  point  are  stored  as  input  quantities.  With  the  information  of  the 
lattice  point  numbers  of  the  tetrahedron  apexes,  the  associated  point 
coordinates  can  be  selected  and  the  stiffness  matrix  [K]  of  the  tetra¬ 
hedron  can  be  computed  according  to  equation  (26).  The  sequence  of  four 
lattice  point  indices  corresponding  to  the  four  apexes  is  computed  and 
stored  for  each  tetrahedron.  In  the  numerical  program,  it  is  called  the 
IV  matrix.  This  matrix  is  later  needed  to  help  identify  the  line  and 
column  of  the  [K]  matrix  with  respect  to  the  displacement  and  fctrce 
indices  (with  the  help  of  equation  47).  After  identification,  a  look-up 
routine  of  the  index  number  in  the  MSK  matrix  allows  generation  of  the 
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matrix  [ K* j ]  appearing  at  the  left-  and  right-hand  sides  of  equation  (37). 
This  look-up  routine  is  a  screening  process  that  is  also  used  to  accumulate 
the  proper  terms  of  the  force  components. 

The  objective  of  this  study  was  to  autogenerate  the  IV  matrix.  This 
was  accomplished  in  two  steps.  The  first  was  to  generate  the  three  lattice 
points  for  each  base  triangle  (the  nonordered  sequence  is  the  IM(J,1) 
matrix  in  the  program,  with  J  -  1,2,3);  the  second  was  to  obtain,  from  the 
nonordered  sequence,  the  ct-@-v  ordered  sequence  of  lattice  numbers  (the 
IM(J,  2)  matrix  in  the  program,  with  J  =>  1,2,3).  The  algorithm  for  both 
steps  follows  from  some  simple  topological  considerations. 

The  configuration  of  triangles  in  Figure  5  allows  ns  to  derive  the 
sequence  of  node  point  indices  as  a  function  of  triangle  index.  Let 
I  5  48  be  the  triangle  index.  Then  we  have  node  1,  the  lattice  index, 


IM(1 , 1)  -  1  + 


[JM 


I-Ie (12)  - 
12 


(48) 


with  Ie (12)  equalling  zero  if  mod(l,12)  =  1  ,  or  one  if  mod(I,12)  =  0  , 


IM(2 , 1)  =  IM(1,1)  +  1  -  (-l)l1/2'j  +  7 


IM(3 , 1)  =  IM(1 , 1)  +  7 


(49) 

(50) 


The  ratios  in  brackets  indicate  that  the  largest  integer  contained  in  the 
quotient  is  taken. 

For  all  triangles  with  numbers  greater  than  48,  a  mirror-image  con¬ 
figuration  exists,  and  the  lattice  indices  are  computed  according  to  the 
following: 


rn  fi-Mi2)' 

L2J  L  12 


IM(1,1)  =  1+|?|+  I - TT— I  (51) 

with  Ie (12)  equalling  zero  if  mod(I,12)  =  1  ,  or  one  if  mod(I,12)  =  0  , 


IM(2, 1) 
IM(3, 1) 


IM(1 , 1)  +  7 


IM(1,1)  + 


(52) 

(53) 


So  far,  we  have  obtained  the  nonordered  sequence  of  lattice  point 
indices  for  each  triangle.  Before  we  continue  to  obtain  the  set  of  a_|3_Y 
ordered  lattice  point  indices,  it  is  appropriate  to  incorporate  the  central 


i 
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symmetry  condition  existing  with  the  problem.  This  symmetry  affects  all 
displacements  at  lattice  points  with  indices  larger  than  32. 

Because  of  central  symmetry  of  the  displacement  distribution,  we  have 


ux(-x,-y,z) 

"  -ux(x,y,z) 

(54) 

uy(-x,-y,z) 

■  -uy(x.y,z) 

(55) 

u2(-x,-y,z) 

=  +uz(x,y,z) 

(56) 

This  means  that  all  displacement  components  at  lattice  points  with  the 
indices  beyond  32  remain  dependent  upon  the  independent  displacement  com¬ 
ponents  already  introduced  at  the  first  32  lattice  points. 

Because  of  equations  (54)  through  (56),  the  axial -symmetric  triangle 
configuration  has  been  chosen  in  Figure  5  and  the  last  considered  lattice 
point  32  terminates  at  the  origin.  The  force  equilibrium  equations  of  all 
components  at  the  first  32  lattice  points  (and  of  the  corresponding  points 
in  the  floors  above)  therefore  represent  one  complete  set  of  independent 
equations  of  the  problem  when  relations  (54)  through  (56)  are  taken  into 
account.  In  order  to  Include  all  [K]  matrix  elements  of  triangles 
associated  with  all  symmetrically  independent  (96  per  floor)  force  com¬ 
ponents,  the  [K]  matrices  up  to  triangle  55  must  be  computed.  Displacement' 
components  for  all  lattice  points  above  32  are  dependent  according  to  relations 
(54)  through  (56).  Therefore,  it  is  not  necessary  to  introduce  more  dis¬ 
placement  indices  or  lattice  point  indices,  but  rather  to  use  the  central 
symmetry  propevty  in  assigning  the  lattice  index  numbers  for  points  beyond 
32.  If,  in  equations  (48)  through  (53),  IM  is  larger  than  32,  it  will 
be  replaced  by  64  -  IM,  and  a  weight  number  IKO  ”  +1  will  be  attached 
to  this  point.  In  this  way,  the  displacement  indices  of  the  symmetric 
points  are  made  equal,  but  the  weight  number  assigned  to  each  lattice 
allows  us  to  distinguish  between  independent  and  symmetrically  dependent 
indices  (and  therefore  displacements).  For  lattice  points  with  index 
numbers  smaller  than  32  (belonging  to  independent  displacements),  the  IKO 
value  is  zero. 


In  this  way,  a  sequence  of  numbers  to  each  base  triangle,  which  is 
called  the  IKO  matrix,  is  generated  along  with  the  IM  matrices.  The 
purpose  of  the  IKO  matrix  is  twofold: 

1.  To  generate  the  coordinates  for  lattice  points  beyond  index  32 
(only  coordinates  of  lattice  points  1  through  32  are  given). 

2.  To  multiply  the  K  matrix  elements  with +1  or  -1  before  collecting 
by  the  summation  process  to  generate  [K*J  .  All  elements 
associated  with  displacement  components  at  lattice  points 
smaller  than  32  are  multiplied  by  +1. 
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[KJmatrix  elements  associated  with  displacements  in  the  z-direction, 
for  lattice  points  beyond  32,  are  also  Multiplied  by  +1  before  collecting. 
However,  the  [Kjmatrix  elements  associated  with  displacements  in  the  x-  and 
y-directions  at  lattice  points  beyond  32  are  to  be  multiplied  by  (-1) 
before  additive  storing  in  the  matrix  (K*l  .  This  scheme  follows  from 
relations  (54)  through  (56). 

For  instance,  in  triangle  44  there  is 


IM(3, 1)  =  64  -  34  -  30 


with  IKO(3,l)  =  1  ;  in  triangle  50  there  is 


IM(3, 1)  -  30 


with  IKO(3, 1)  -  0  . 

With  the  IM(J , 1)  matrix  and  the  IKO(J,l)  matrix  generated,  where 
J  ■  1,2,3  ,  the  lattice  points  for  a  specified  triangle  are  known.  Since 
the  index  number  of  the  lattice  point  in  the  IM  matrix  does  not  reveal 
whether  it  is  an  independent  or  axial-symmetric  dependent  point,  the  IKO 
matrix  carries  this  information  for  this  purpose. 


From  the  IM(J,1)  matrix,  we  get  the  a-0-Y  ordered  matrix  by  using 
the  topological  pattern  existing  for  the  lattice  points  in  each  row.  The 
first  row  points  alternate  between  y-a  points,  the  second  row  points 
between  P  -y  points,  and  the  third  row  points  between  a-P  points;  the 
fourth  row  shows  the  same  cycling  as  the  first  row,  and  the  fifth  shows  the 
same  cycling  as  the  second  row.  We  assign  to  each  point  with  the  lattice 
number  I  a  number  ©(I)  ,  such  that 


0(1) 


0  (i) 


0(D 


!if  I  is  a 
topological 
a-point 


{if  I  is  a 
topological 
0-point 


{if  I  is  a 
topological 
y-point 


Corresponding  to  the  y-a  alternation  of  the  first  seven  lattice 
points,  the  associated  0  values  form  a  sequence  0, +1,0, +  1,0, +  1,0  .  The 
next  seven  form  a  sequence  -1, 0, -1, 0, -1, 0, -1  ,  and  the  seven  lattice 
points  in  the  third  rc w  form  a  sequence  +1 , -1 ,+l , -1  ,+l , -1  ,+l  . 
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The  following  21  numbers  are  periodic  to  the  first  21.  This  sequence 
of  21  numbers  is  identical  with  the  elements  of  the  sequence  {  0  (I)  }  , 
where  I  ■  1,2,., 21  ,  with 

0  (I)  =  3  {  1+'^1  [  2  cos  *  (1  +  2k)  -  2  cos  3  (3  +  2k)] 

-  (-Dk  [l  -  2  cos  3  (-1  +  2k)]  |  (57) 

with  k  =  [  1-Ie (7) / 7 ]  and  with  Ie(7)  equalling  one  if  mod(I,7)  =  0  , 
or  equalling  zero  if  mod(I,7)  i-  0  .  Formula  (57)  serv'es  as  a  guide  in 
making  a  decision  about  the  topological  property  of  the  lattice  point 
having  the  index  I  .  Figure  7  is  a  flow  diagram  of  a  version  following 
the  logical  content  of  formula  (57)  that  was  used  in  this  investigation. 


Figure  7.  Computer  Flow  Diagram  to  Determine  Topological 
Property  of  Lattice  Point. 
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Going  with  each  lattice  point  index  through  this  routine,  we  generate 
the  a-0'Y  ordered  vector  {IM(J,3}  ,  where  J  =  1,2,  3  ,  with  the  meaning 


IM (1,2)  = 

index  of 

a-poinL  ^ 

of  a 

IM(2, 2)  = 

index  of 

B-point  > 

particular 

base 

(58) 

IM(3,2)  = 

index  of 

y-poin t  ) 

triangle 

With  the  {1M|  matrix. 

the  four 

lattice  point  indices  of  the 

apexes 

all  three  tetrahedrons  above 

the  base 

triangle 

follow  immediately 

from  the 

above  definition  of  generating  the  apexes  of  the  tetrahedrons. 

First  tetrahedron  above 

the 

triangle 

: 

IV  (J, 

1,0  = 

IM(1,2) 

IV(J, 

2,D  = 

IM(2 ,2) 

IV(J, 

3, 1)  = 

IM(3, 2) 

IV  (J, 

4,1)  = 

IM(2 , 2)  + 

32 

(59) 

Second  triangle  above  the  triangle: 


IV(J , 1 , 2) 
IV(J,2,2) 
IV(J, 3,2) 
IV(J,4,2) 


IM(1,2) 
1M(2,2)  +  32 
IM(3, 2)  r  32 
IM(3, 2) 


th 


Third  tetrahedron  above  the  J  triangle: 


(60) 


IV(J , 1 , 3) 
IV  (J ,  2 , 3) 
IV(J, 3, 3) 
IV (J ,4 , 3) 


IM(1,2)  +  32 
IM(2,2)  +  32 
IM(3 ,2)  +  32 
IM(1, 2) 


(61) 
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With  the*  IV  matrix,  the  set  of  four  apex  indices  is  known;  therefore, 
the  coordinates  of  the  apexes  can  he  obtained  and  the  elements  of  the  [K] 
matrix  by  means  of  equation  (30)  can  now  be  computed.  The  IV  matrix  must 
be  saved,  since  it  is  used  to  identify  row  and  column  of  the  [k]  matrix 
as  well  as  to  generate  the  desired  elements  of  the  [ K J  matrix. 

In  the  next  section,  a  description  is  given  of  how  the  desired 
elements  of  the  matrix  [K*]  are  generated  from  the  matrices  [  K  J  , 


Generation  of  the  Matrices  Required  To  Solve  for 
the  Unknown  Displacements 

The  main  objective  was  to  find  the  coefficient  matrices  of  equation 
(37)  and  to  find  solutions  for  the  unkown  displacement  components  u'u'  , 

where  j  -  M  +  1  ,  M  +  2,  . .  .M.  J 

l  l 

Earlier  in  this  report,  a  detailed  description  was  given  of  the  gener¬ 
ation  of  the  elements  of  the  matrix  £ K*J  .  Equation  (31)  reveals  the 
algorithm  for  obtaining  the  coefficients  of  £ K*J  ,  and  equation  (37) 
indicates  which  of  the  coefficients  are  actually  used.  These  two  prop¬ 
ositions  are  observed  in  generating  the  matrix  elements  of  equation  (37). 

The  actual  procedure  calls  for  re-indexing  of  the  unknown  displace¬ 
ments  in  a  consecutive,  natural  number  sequence,  since  the  indexing  derived 
from  the  lattice  point  indices  (see  equation  47)  represents  a  nonsequen- 
tially  ordered  subset  of  natural  numbers.  Re-indexing  is  conveniently 
done  by  means  of  the  sequence  containing  the  indices  of  unknown  displace¬ 
ment  components  (or  known  force  components)  as  elements.  This  sequence  is 
called  the  MSK  matrix  and  has  been  generated  simply  by  picking  up  only 
displacement  indices  which  are  unknown  in  accordance  with  boundary  condi¬ 
tions  (38)  through  (45).  The  first  three  numbers  in  the  first  column  at 
each  lattice  point  are  the  indices  of  the  x,y,z  displacements  at  that  lattice 
point.  Behind  each  index  is  the  argument  number  of  its  appearance  in  the 
MSK  matrix. 

From  the  three-dimensional  lattice  configuration,  wherein  Figure  5 
is  the  ground-floor  portion,  we  can  see  that  the  first  80  equilibrium 
equations  contain  terras  involving  the  first  160  unknowns  only.  This  is 
because  the  force  components  of  the  first  floor  are  influenced  by  dis¬ 
placement  coefficients  of  the  first  and  second  floors  only.  The  first 
30  lines  of  the  matrix  (K*j  thus  have  nonzero  elements  in  the  first 
160  columns  only;  all  further  columns  contain  zero  elements. 

The  80  equilibrium  equations  of  all  known  force  components  at  the 
second  floor  (equations  81  through  160)  are  in  general  connected  to  un¬ 
known  displacements  in  the  first,  second,  and  third  floors.  This  means 
that  the  [K*]  matrix  elements  for  lines  81  through  160  contain,  in  general, 
nonzero  elements  in  the  first  240  columns. 
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The  equilibrium  equations  of  the  third-floor  force  components  are  of 
the  same  nature  as  the  corresponding  equations  of  the  second  floor.  The 
only  difference  is  that  all  the  unknown  displacement  coefficient  indices 
increase  by  80.  This  means  that  the  submatrix  contained  in  lines  81 
through  160  is  repeated  in  lines  161  through  240,  but  that  all  columns  are 
shifted  to  the  right  by  80  columns.  The  same  considerations  apply  for  all 
floors  following.  Matrix  [K*]  therefore  contains  three  diagonal  non¬ 
zero  submatrices  which  are  the  same  in  the  entire  matrix,  except  for  the 
first  and  last  diagonal  matrices.  Figure  8  depicts  the  [ K*]  matrix. 

Because  the  [K*]  matrix  is  symmetric,  only  three  different  types 
of  80  X  80  submatrices  are  involved  in  defining  [K*]  .  These  are 
the  submatrices  [A],[b],(C]  .  This  means  that  it  is  sufficient  to  gen¬ 
erate  the  [K*]  elements  for  the  first  160  lines  only.  In  other  words, 
only  the  proper  elements  of  the  f K]  matrices  of  the  set  of  tetrahedrons 

between  the  first  and  third  floors  need  to  be  taken  into  account  in  order  to 
generate  [A),[b],[C]  . 

The  scheme  applied  to  generate  the  matrix  elements  of  the  80  X  80 
matrix  [AS22]  -  j_[A],[B]^  and  the  80  X  160  matrix  [KS22]  -  Q[B ] , 
[CJ,[B]]  is  basically  a  screening  routine  applied  to  each  [K J  matrix 
element  extended  over  all  [K]  matrices  of  tetrahedrons  in  the  first  two 
layers.  The  [K]  matrix  of  a  particular  tetrahedron  in  the  first  layer  is 
identical  to  the  [K]  matrix 
of  the  tetrahedron  that  is 
its  second-layer  counterpart. 

Therefore,  the  ( K]  matrices 
for  the  three  first-floor 
tetrahedrons  above  a  tri¬ 
angle  are suf f icient,  for  our 
purposes,  to  generate  [AS22] 
and  t  KS  22]  . 


In  the  process  of  giving 
a  program  description,  it  was 
shown  in  the  preceding  sec¬ 
tion  that  the  [ K J  matrix 
for  each  of  the  three  tetra¬ 
hedrons  above  a  base  triangle 
is  obtained.  With  the  addi¬ 
tional  information  gathered 
in  this  section,  it  is  now 
possible  to  identify  each 
element  of  the  matrix  with 
respect  to  that  element  of 
the  matrix  [AS22]  or  [ KS 22] 
(including  the  right-hand - 
side  independent  terms)  to 
which  it  will  be  added.  A 
diagram  of  this  scheme,  Fig¬ 
ure  9,  shows  how  this  is  done. 
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Figure  9.  Computer  Flow  Diagram  for  Calculating  The  Stiffness  Matrix  [K*J . 
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Since  the  IV  matrix  is  generated  for  each  triangle  index  I  according 
to  Figure  9,  the  oniy  input  required  for  this  routine  is  the  x-  and  y- 
coordinates  of  the  32  ground-floor  lattice  points. 

After  compl*  ion  of  all  the  indicated  looping  processes,  all  of  the 
12  X  12  elements  of  each  [K]  matrix  belonging  to  every  tetrahedron 
in  the  two  layers  have  been  screened  and  properly  accumulated  in  the 
matrices  [AS22]  and  [ KS22 ]  as  well  as  in  the  independent  vectors 
[A21]  and  [K21]  .  Therefore,  the  elements  of  the  matrices  [AS22]  and 
[KS22]  along  with  the  independent  vectors  [ A2 1 ]  and  [K21]  are  complete 
and  considered  to  be  generated.  The  matrix  represented  in  Figure  8  is 
therefore  established. 

The  systems  of  equations  belonging  to  the  Figure  8  matrix  is  undeter¬ 
mined,  since  the  matrix  has  720  lines  and  800  columns.  The  80  missing 
equations  are  obtained  by  taking  the  internal  boundary  conditions  (43), 
(44),  and  (45)  —  in  other  words,  those  "inside"  the  body  —  into  account. 
In  line  with  the  first-order  approximation  used  throughout  this  study,  we 
can  enforce  an  approximate  version  of  (43)  and  (44)  by  making  all  u^  and 

Uy  components  at  the  tenth  floor  equal  to  the  corresponding  u^  and  u^ 

components  at  the  ninth  floor.  By  corresponding  u  and  u  components,  we 

x  y 

mean  that  the  u  and  u  values  for  lattice  points  with  the  same  x-  and 
x  y 

y-coordinates  in  both  floors  are  considered.  Condition  (45)  is  satisfied 
by  putting  all  u z  values  equal  to  zero  at  floor  number  10.  These  three 

additional  conditions  represent  80  linear  equations  which  now  supplement 
the  set  of  720  linear  equations  with  800  unknowns  (reference  the  Figure  8 
matrix).  The  80  equations  are  as  follows: 

Un(i)+720 
Um(j)+720 

where  i  =  1,2... 48  and  j  =  1,2... 32  .  Here,  n(i)  is  the  sub-sequence 
of  the  first  80  elements  of  the [MSK]  matr ix  that  contains  only  x-  and  y- 
directed  unknown  displacement  indices;  m(j)  is  the  sub-sequence  that 
contains  elements  of  the  [MSK]  matr ix  associated  with  unknown  z-directed 
displacement  indices. 

The  unknown  with  indices  greater  than  720  appear  in  the  last  80 
equations  of  the  system  of  720  equations.  Matrix  [B]  represents  the 
coefficients  for  them.  Substitution  of  equations  (62)  and  (63)  into 
these  last  80  equations  modifies  the  set  of  equations  to  the  extent  that 
matrix  [B]  is  absorbed  into  matrix  [ C J  . 

We  may  write  the  last  80  equations  in  the  following  manner. 


n(i)+640 


(62) 

(63) 
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u(u)+ 
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48 
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K*  >> 

i,640+n(j)  640+  n(j) 


32 

+V*  K* 

i  ,  640+  m  ( j ) 


U64u+m  (j  ) 


48 

+V'  K*  (u) 

4“f  i ,  720+n  ( j  )  U720+n(j) 


32 

+V  K*  u(u) 

i,720+m(j)  720+m  (j) 


(64) 


where  i  =  641,  642,  ...  720  .  In  equation  (63) , "we  have  the  first  sum 
term,  the  [B]  matrix  term  of  Figure  8.  The  next  two  terms  (the  [C] 
matrix  terms),  the  last  two  terms  (the  [Bj  matrix  terms),  and  the  two 
separate  terms  for  the  [B  ]  and  [C  ]  matrices  correspond  to  the  separa¬ 
tion  into  x-  and  y-directed  displacements  and  z-directed  displacements. 
Now  we  substitute  equations  (61)  and  (62)  into  (63)  and  obtain  an  equa¬ 
tion  containing  only  displacement  indices  561. ..720  : 


640 

Y  K*.  u<u) 

iitii  3  3 


48 

^  (Ki,640fn(j)  +  Ki,720fn(j)/ 
32 

+Y  K*  u(u) 

4^  i,640fm(j)  640fm(j) 


u 


(u) 

64  Of  n(j) 


(65) 
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where  i  =  641,  642,  ...  720  .  The  first  term  (corresponding  to  the  [B] 
matrix  and  the  independent  term)  is  not  altered.  The  second  two  sum  terms 
correspond  to  the  new  matrix  [ C ' J  . 

According  to  equation  (65),  matrix  [ C ' ]  can  be  generated  in  the 
following  manner.  First,  all  column  elements  of  the  [b]  matrix  with  the 
indices  m(j)  ,  where  j  =  1,2,... 32  ,  are  zeroed  out.  The  matrix  obtained 
is  called  [B1]  and  is  an  (80  x  80)  matrix.  Matrix  (C'J  is  the  sum  of 
matrices  [C]  and  [B'J  . 

With  the  "internal"  boundary  conditions  (61)  and  (62),  along  with  the 
system  of  equations  based  on  the  matrix  depicted  earlier  in  Figure  8,  we 
can  obtain  an  equivalent  system  or  720  equations  which  correspond  to 
the  matrix  inside  the  broken  line  of  Figure  8.  In  this  case,  the  sub¬ 
matrix  [C]  of  the  last  line  is  replaced  by  the  matrix  [C'J 


Solutions  for  this  system  of  equations  have  been  obtained  by  apply¬ 
ing  the  method  of  Choleski.*  The  accuracy  of  this  method  has  been 
verified  by  computing  a  test  case  for  homogeneous  material,  so  that 
the  stress  and  displacement  distribution  could  be  obtained  by  analytical 
means . 

For  our  boundary  problem,  the  displacements  in  the  homogeneous  test 
case  are  exactly  linear.  Therefore,  the  approximation  used  in  this  program 
becomes  an  exact  solution.  The  displacements  obtained  as  solutions  of 
the  system  of  720  linear  equations  (using  the  Choleski  method)  agreed 
in  the  first  7  digits  with  the  displacements  obtained  by  analytical 
means . 

This  test  result  carries  with  it  the  implication  of  the  correctness 
of  the  IK*]  matrix  coefficients.  Correctness  of  the  [K*]  matrix  is 
maintained  for  the  fiber-matrix  case,  since  the  scheme  to  generate 
the  elements  of  [K*]  is  independent  of  input  values. 


Computation  of  the  Stress  Components 

The  stresses  for  each  tetrahedron  can  now  be  calculated  with  the 
use  of  equation  (21),  since  all  u's  are  now  known;  in  other  words, 
they  are  either  given  or  are  solutions  of  the  u^u)  equations. 

The  remaining  task  is  to  find  the  12  proper  displacement  components 
in  the  correct  order  for  each  tetrahedron  in  every  layer.  This  is 
accomplished  in  a  manner  similar  to  the  way  the  equation  indices  of  a 
Ik]  matrix  element  were  obtained. 


*  Also  known  as  Crout's  scheme. 


r 
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In  Figure  10,  the  scheme  to  obtain  the  six  stress  elements  (the 
a  matrix)  for  each  tetrahedron  is  demonstrated.  Given  are  all  known 
displacements  and  the  solution  vector  u(i)  ,  where  i  =  1,2,... 720  . 

The  meaning  of  the  six  stress  components  in  the  order  appearing  in 
equation  (16)  is  obtainable  from  equation  (10): 
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With  all  of  the  a  components  obtained,  the  stress  distribution  within 
the  three-dimensional  domain  can  be  obtained, since  the  coordinates  of  the 
tetrahedron  centroids  are  the  field  point  coordinates  of  this  distribution 
approximation. 
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Figure  10.  Computer  Flow  Diagram  for  Calculating  Stress  Components. 
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Select  estimated  argument  of  MSK.  matrix  such  that 
associated  MSK  element  is  as  close  as  possible  to 
the  displacement'  index  I,j 


If  estimation  was  erroneous,  iterate  on  argument  of 
MSK  matrix  until  difference  between  MSK  element  and 
lcl  is  zero  or  shows  change  of  sign.  Zero  means 
is  e'ement  of  MSK  matrix;  sign  change  means  Ij  is 
not  element  of  MSK  matrix  (unknown  or  known  dis¬ 
placement  )  . 


is  contai ni 
in  MS K  matrix 


gument  of  MSK  matrix  such  that 
nent  is  as  close  as  possible  to 
idex  Ij 


•roneous,  iterate  on  argument  of 
:ference  between  MSK  clement  and 
change  of  sign.  Zero  means  Ij 
itrix;  sign  change  means  I<]  is 
latrix  (unknown  or  known  dis- 


is  not  contained 
in  MSK  matrix 


After  exhaustion  of  all  Kj  and  Ka , 
all  12  displacements  f>x(Kz),  wht* i*t* 
K  =1,2,...  12,  of  tetrahedron  I3  of 
triangle  Ij  at  Floor  1  are  obtained 
in  proper  order. 

Now  we  can  calculate  the  stresses 
for  this  particular  tetrahedron 
according  to  equation  (16),  since 
the  matrices  ( C  ] [ D ] |  A ] _1  =  {CDA} 
have  been  obtained  as  a  by-product 
when  matrices  {k)  were  calculated 
and  then  stored. 


a  =  |cda| J  dx 


The  six  a  values  are  printed  along 
with  the  coordinates  of  the  centroid 
of  the  tetrahedron 


NUMERICAL  RESULTS  FOR  THE  THREE-DIMENSIONAL  SOLUTION 


Figure  11  depicts  a  cross  section  of  a  composite  under  transverse 
loading.  Figures  12  through  20  illustrate  the  variation  of  the  stresses 
in  the  fiber  direction  at  the  interface  points,  indicated  in  Figure  li, 
in  regions  near  the  free  end  of  the  composite.  These  stresses,  which 
appear  because  of  the  difference  of  elastic  constants  in  the  fiber  and 
the  resin,  are  similar  to  the  stresses  obtained  with  the  plane  stress 
assumption  for  planes  far  from  the  free  end. 

The  three-dimensional  analysis  shows  the  existence  of  shear  stresses 
CTxz  anc*  °\z  considerable  magnitude  at  the  interface  points  depicted 

in  Figure  11.  In  fact,  these  shear  stresses  are  about  one-half  the  peak 
shear  stress  aXy  from  the  plane  analysis.  The  numbers  have  been  derived 
from  numerical  calculations  of  the  stresses  within  the  area  of  triangle 
41  (see  Figure  5) . 


Figure  11.  Cross  Section  of  a  Composite 
Under  Transverse  Loading. 
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Figure  12.  Longitudinal  Stresses 
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Figure  13.  Longitudinal  Stresses 
ax  at  L|ocation  Indic¬ 
ated  in  Figure  11 
Due  to  Transverse 
Normal  Loading  ax 
Average . 
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Figure  14.  Longitudinal  Shear 
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Location  Indicated 
in  Figure  11  Due 
to  Transverse  Norm¬ 
al  Loading  ax 
Average. 
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Figure  18.  Longitudinal  Shear 
Stresses  aXz  at 
jLocation  Indicated 
in  Figure  II  Due 
to  Transverse  Normal 
Loading  Average. 
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Longitudinal  Shear 
Stresses  oxz  at 
Location  Indicated 
in  Figure  11  Due 
to  Transverse  Normal 
Loading  Average. 


Figure  21. 
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SOLUTIONS  OF  THE  PLANE  PROBLEMS 


In  regions  of  the  composite  under  transverse  loading  that  are  far 
from  the  ends  (see  Figure  22),  it  is  possible  to  realistically  assume 
plane  strain  or  plane  stress.  It  will  be  plane  strain  if  the  displace¬ 
ment  u2  is  constrained  at  the  ends  z  =  0  and  z  =  ;  it  will  be 

plane  stress  if  the  ends  are  free  of  loads. 


Figure  22.  Composite  Under 
Transverse  Load. 


It  is  assumed  that  the  transverse  loads  P  are  applied  to  the  com¬ 
posite  far  from  the  chosen  representative  element  of  analysis.  Also,  it 
is  assumed  that  the  composite  cannot  have  expansions  in  the  y-direction; 
in  other  words,  compared  with  the  dimensions  of  the  typical  element,  the 
length  of  the  composite  is  infinite  in  the  y-direction. 
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With  these  assumptions, 
element  (Figure  23)  are: 


the  boundary  conditions  for  the 


represents  tive 


Figure  23.  Representative  Element. 


At  the  interface,  between  the  fiber  and  matrix,  the  following  conti¬ 
nuity  conditions  must  be  satisfied: 
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nt 


(67) 


The  superscript  f  indicates  fiber  and  the  superscript  m  indicates 
matrix.  The  n  and  t  are  the  normal  and  tangential  directions  at  the 
interface  points,  respectively. 

Moreover,  a  symmetry  condition  exists  for  the  displacements  and,  con¬ 
sequently,  for  the  stresses  of  the  representative  element.  That  is,' 


u  (x,y)  =  -u  (-x,-y) 

X  X 

uy(x,y)  =  -uy(-x,-y) 


(68) 


The  stress  distribution  was  solved  using  the  finite  element  method 
with  the  triangular  net  shown  in  Figure  24.  This  was  explained  in  detail 
for  the  three-dimensional  solutions  in  the  last  section. 
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Figure  24.  Finite  Elements  of  a  Basic  Representative  Cross-Sectional  Area 
of  Floor  Number  1  Used  in  Two-Dimensional  Numerical  Analysis. 


Once  the  stress  distribution  is  known,  the  transverse  modulus  of  the 
composite  can  be  found.  In  fact, 


where 


Then, 


c  -  fi 


b 

2 


P 

2  '  k 


(69) 


The  total  applied  force  P  is  obtained  by  adding  the  forces  at  the 
x  direction  (ctx  •  dimension  at  the  y  direction)  for  the  triangles 
bordering  the  boundary  x  =  c  . 

By  divising  equation  (69)  setting  k  =  1  by  the  modulus  of  the  fiber 
Ej  ,  the  nondimen sional  expression  is  obtained. 


X 


(70) 


In  the  following  section,  numerical  results  of  the  two-dimensional 
analysis  are  given. 
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NUMERICAL  RESULTS  OF  THE  PLANE  PROBLEMS  AND  TEST  RESULTS 


The  trajectories  of  the  principal  stresses  in  the  typical  element  of 
a  composite  material  under  transverse  loading  are  represented  in  Figure  25. 
The  nonhomogeneity  introduced  by  the  presence  of  the  fiber  causes  a  devia¬ 
tion  in  the  trajectories  with  respect  to  the  homogeneous  case  where  the 
trajectories  are  obviously  parallel  to  the  x  and  y  axes.  This  fact 
suggests  the  presence  of  stress  concentrations,  which  will  be  demonstrated 
effectively  in  the  following  discussion. 

A  parametric  study  was  performed  to  evaluate  the  influence  of  the 
matrix  and  fiber  properties  on  the  stress  distribution.  Figure  26  shows 
the  radial  and  tangential  stresses  along  the  interface  for  composites  with 
60  percent  of  the  fibers  containing  the  same  matrix  material  but  with 
different  fiber  material.  One  composite  has  a  modulus  relationship  of 
Ef/Em  =  20,  and  the  other  has  a  relationship  of  Ef /Em  =  120  .  These 
correspond  approximately  to  glass  fiber  and  boron  fiber  composites  with 
epoxy  resin,  respectively.  From  the  curves  of  the  figure,  it  is  possible 
to  deduct  the  slight  influence  the  material  of  the  fiber  has  in  the  stress 
distribution.  This  conclusion  is  correct  in  all  cases  In  which  the  fiber 
is  considerably  harder  in  comparison  with  the  matrix,  as  usually  occurs  in 
most  of  the  composites. 

The  nondimensional  ordinates  o_/a  and  T  n/<7  of  the  curves 

r  avg  ry  avg 

were  obtained  by  dividing  the  actual  stresses  ar  ana  Trg  by  the  average 
stress  caVg  found  by  averaging  the  stresses  <jx  at  the  triangles  3,  12, 
24,  36,  44  and  49. 

Another  parametric  study  was  performed  by  taking  the  volumetric  con¬ 
tent  as  variable,  keeping  tne  modulus  relationship  Ef/Em  constant.  The 

results  are  indicated  in  Figures  27  and  28,  where  the  stresses  along  the 
interface  are  plotted.  As  is  easy  to  imagine,  the  peak  stresses  are  greater 
when  Vp  is  increased. 

Figure  29  indicates  the  displacement  component  at  the  x-direction  for 
composites  with  Ef  /Em  ■  20  and  120  .  In  these  curves  it  is  possible  to 
observe  the  small  influence  of  the  fiber  deformation  in  comparison  with  the 
total  deformation,  which  is  carried  almost  completely  by  the  matrix. 

TRANSVERSE  MODULUS 


A  parametric  study  was  performed  to  compute  the  transverse  modulus  of 
several  composite  types.  Two  different  conditions  were  considered:  (1) 
plane  stress  on  the  x,y  plane,  and  (2)  plane  strain  on  the  same  plane. 


45 


.« 


Figure  25.  Stress  Trajectories  of  Transverse  Loading. 
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Figure  27.  Stresses  Along  the  Interface  With 
Volumetric  Content  Variable  and 
Modulus  Relationship  Constant. 
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Figure  29.  Displacement  Component  at  the  x-Direction. 
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Figure  30  gives  the  transverse  modulus  of  the  composite  E^.  for  a 
plane  stress  condition.  Figure  31  gives  the  modulus  of  the  composite  for 
plane  strain  as  a  function  of  the  modulus  ratio  ,  and  by  taking  the 

volumetric  content  as  parameter.  It  can  be  appreciated  that  the  influence 
of  the  fiber  modulus  on  the  composite  modulus  is  small,  in  accordance  with 
the  displacement  distribution  discussed  previously. 

Figure  32  gives  a  comparison  of  the  transverse  modulus  obtained  with 
Ekvall' s  formula, 


and  the  results  from  the  computer  analysis. 


m 

Figure.  30.  Transverse  Modulus  of  a  Composite  as  a 


Function  of  Fiber  and  Matrix  Modulus  and 
Volume  Percentage. 
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Figure  33  presents  the  comparison  with  reppect  to  tht  Shaffer's 
formula, 


(72) 


and 


for  Vf^0.68  (73) 


As  observed  in  the  figures,  the  transverse  modulus  from  the  micro¬ 
mechanical  analysis  is  greater  than  the  corresponding  modulus  in  the 
Ekvall  and  Shaffer  formulas,  especially  for  the  higher  volumetric  contents. 

The  transverse  modulus  of  a  composite  obtained  by  plane  stress  analy¬ 
sis  was  presented  in  Figure  30  for  different  volumes  and  material  contents 
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of  the  reinforcement  and  the  matrix.  This  modulus  can  be  approximated 
by  the  following  formula  : 


E 

in 


(74) 


It  is  interesting  to  note  that  this  formula  is  identical  to  Rosen's 

formula*  for  longitudinal  shear.  The  difference  is  that  instead  of  the 

G  moduli,  the  E  moduli  of  the  different  materials  must  be  used  as  indicated 

in  equation  (74).  Allowing  for  little  error  for  low  fiber  percentages 

and  E./E  <100  ,  this  formula  has  an  error  showing  a  4  percent  higher  trans- 
f  m 

verse  modulus  at  higher  values.  The  following  is  a  formula  which  permits 
errors  of  less  than  1  percent  for  all  possible  combinations. 
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6*666Vf 


1-27 


(75) 


*  B.  W.  Rosen,  N.  F.  Don,  and  Z.  Hashin,  Mechanical  Properties  of  Fibrous 
Composites ,  NASA  CR-31,  Contract  NAS-470,  General  Electric  Corporation, 
April  1964. 
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APPENDIX  I 

PROBLEM  OF  FIBER- REINFORCED  COMPOSITE  SUBJECTED  TO 
TRANSVERSE  LOADING  SOLVED  BY  POINT- MATCHING  METHOD 


When  the  weight  is  the  only  body  force  ,  the  plane  strain  problem  can 
be  solved  by  finding  the  polynomial  solution  of  the  differential  equation 


£_£  +  2  — &lfi_  +  Sil  =  o  (76) 

ax4  ay4 


which  satisfies  the  boundary  conditions.  In  equation  (76),  0  is  the 

Airy  stress  function.  Stresses  and  displacements  are  then  defined  as 
follows  : 
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.  vM] 

ay  J 

+  f2(x) 

(81) 

The  strains  are 


du 

w  x 

ax 

au 

__x 

ay 


l+y 

E 

1+v 

E 


H  &  ■  ’  0] 

[. v  ail  +  (1_v)  sill 
L  ay8  a*2  J 


"xy 


au  au 
— + 
ay  ax 


zq+y)  as 

e  axay 


(82) 

(83) 

(84) 
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From  equations  (80),  (81),  and  (84),  we  have 


|dx+  r^dy+2^+-i- 

J  hy  J  ax3  sxay  l-va 


rdf.  df. 


v2  l  dy  dx 


+  -rrl  -  o  (85) 


The  polynomial  solution  to  the  Laplace  equation  (harmonic  equation)  is 


(x+iy) 


I 

n'  =  0 
N 

I 

n=0 

N 

I 

n=0 


N  1  N n*  N-n' 

(iy)  x 


. 2n  N-2n 
(iy)  x  + 


E(> 


..  .2rrhl  N-2n- 1 
(iy)  x 


N 

,  ,Nn  2n  N-2n  .  V  I  N  \  ,  ,.n  2n+l  N-2n-l 
(-1)  y  x  +  i  )  (-1)  y  x 

n=0 


>N  (x,y)  +  i  xN(x,y) 


(86) 


Therefore,  the  solution  to  the  biharmonic  equation,  in  region  f (fiber)  , 
is  as  follows: 


6  =  K^xy9  +  Kxya  +  K3xy  +  K^x  +  KBy3  + 


Kay3  +  V  +  V3y  +  K9x3y  +  K10x3  +  K1Tx  = 


M  m 

2ml  /  i \n  2n  2m-2rH-l 
,  ,  (-1)  y  x  + 

m  4-  \2n/ 

10=2  n=0 


I  I 


M 


m 


.  ,n  2rrfl  2m-2u 
B  )  (-1)  y  x  + 

»  1 2n  ,* 

trf=2  n=0 


I  ■-  I 


M  m 

i  t  :;:i 

_ n  _ n  • 


.  . . n  2n+l  2m-2n+l  , 
(-1)  y  x  + 


mp2  n=0 
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t 


M  m  , 

£■■£  O 


,  ,.n  2n+2  2:n-2n 

(-1)  y  x 


npl  n“0 


By  substituting  g  from  equation  (87)  into  equations  (77)  through  (81), 
the  following  stresses  and  displacements  are  obtained: 


6KTxy  +  2Kgx  +  6KBy  +  2K  + 


M  m 

Y.  *•  L 0  <2n)(2-i>(-i,n  *2<" 

_ n  _ a  '  ' 


-1)  2m-2trfl  . 
x  + 


EC) 


(2n+l)(2n)(-l)n  y2n-1  x2m-2n  + 


w  ,\n  2n-l  2m-2n+l 
(2n+l)(2n)(-l)  y  x  + 


£(;;:) 


(2n+2)(2n+l)(-l)n  y2n  x2m'2n 


6K8xy  +  2K#y  +  6K1Qx  +  2KJX  + 


eCi) 

n°0  ' 

to 


(2m-2n+l)(2m-2n)(-l)n  y2"  x2m'2n_1  + 


(2m-2n)(2m-2n-l)(-l)n  y2n+1  x2m'm"2  + 


m 

^  (2m-2n+l)(2m-2n)(-l)n  y2n+l  x2m-2n-l  + 


(Continued) 


it 


(89) 


t  D"  Ifcj  <2"-2n)(2"- 

nrl  n-0  ' 


2n-l)(-l)n  y2n+2  x2m'2n-2 


'xy 


-3KlYa  -  2Kay  -  K3  -  3KflxS  -  2K9x  - 
M  m 


I  V  l 

nr  2  n“0 

M  m 


I  B-  I 


l2rW 
f  2m\ 


(2n)(2m-2n+l)(-l)n  y2n_1  x2m'2n 


B„  ^  II  (2n+l)(2m-2n)(-l)n  y2n  x2m_2n'1 


nr  2  n=0 

M  m 

C 

nr  2  n*0 

M  m 


M  m 

y  cm  y  y2n  x2m-2n . 

_  r>  /»  *  ■ 


Y  Dm  Y*  2I+I  (2n+2) (2m-2n) (-l)n  y2n+1  x2"1'2"*1 

1  __n 


f 

“x 


nr  1  n“0 


1+v  ( 

—  -  ^[3  (l-vfjxay  -  (2-vf)y3]+  K3[(l-vf]x®  - 
(2-vf)y»]  -  K3(2-vfjy  -  K4vf  +  6KB(l-vf)Xy  + 

a.hf|x  -  Kb [j1_vf)y3  +  3Vay]  -  (2vf^y)  - 

K10[3  J-Vf)ya  +  3vfx8]  -  +  w0y(1-vf]}  + 

1  I  ^1  M  rn 

T  [  Q  <“>"  [M1!^  • 

(  nr  2  n=  0  L 

y2(n-l)  x2m-2rH-2  _  Vf(2m_2irfl)  y2n  y2m-2nj  + 


(90) 
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M  m 

+  V  Bm 

m“2  n=0 


z  b»  z  (;j  <-»■  m 

n=0 

2n-l  2m-2n+l  /0  „  .  2n+l  2m-2r.-l1 

x  -  v^(2m-2n)  y  x  j  + 


11  m 

Zr— <  2nrl-l 

Z| 


m=2  n=0 

2n-l  2m-2rr+2 

y  x 

M 


2n+lJ 


1  '  [\  f|  2m-2n+2 


1  . 


/">  o~li  s  2n+1  2m 
-  Vf(2m-2n+l)  y  x 


-2„]  + 


n  in 

y  D  y  2"*1  (.i)n  rii.v.i  . 

L-,  m  Zj  2n+l|  1  [l  fl  2m-2n+l 


nr  1  n=0 


L+v 


2n  2m-2n+l  ,  .  2n+2  2m-2n- 

y  x  -  vf(2m-2n)  y  x 


{ -Ki[[ivf)x3  +  -  |ta|2v*xy)  ■  k>M' 

[k,  3|l-vf]x>+  3-Jfy=]  -  Ks(2v£y]  -  K,vf  + 

K,[3(l-Vf)xy>  -  (2-V£)x=]  +  K,[(l-vf]y*  -  (2-vf)x»]  + 
[jl-vf)6xy]  +  Ku[2(l-Vf)y]  +  W„[- [l-.f)x]J+ 

Mm. 

z  z  3  <-i>“  M  i2m~2^<2*~2n) 


10 


l+v 


nr2  n=0 


2rr+l  2m-2n-l  .  2n-l  2m-2n+ll 

y  x  -  Vf(2n)  y  x  J  + 


M 


m  / 

z  b»  z  d  <-i>n  im  (2"r2^2m~2"~i) 


nr  2  n=0 


J 


(Continued) 
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I 


2n+2  2m-2ri-2  .  2n  2m-2nl 

y  x  -  vf(2n+l)  y  x  J  + 


M  m  .  . 

I «.  iD  <-»■  h  ^--£&"~2ni 


m=2  n“0 


2n+2  2m-2n-l  2n  2m-2rr(-ll 

y  x  -  v^(2n+l)  y  x  J  + 


M  m 

E  I  □  (-1)n  M 


np  1  rv=0 


2n+3  2m-2n-2  .  2n+l  2m-2nl 

X  -  V^(2n+2)  y  x 


(92) 


In  order  to  satisfy  the  equation  fi(x,y)  =  Q(-x,-y)  ,  the  polynomial 
solution  to  the  biharmonic  equation  in  region  m(matrix)  is  simplified 
as  follows : 


0m  -  Kxy3 


M 


r  3 


z 

mf=2 


,xy  + 

Va  + 

Kax  y  +  Ku 

m 

y 

/  2irr+- 1  \ 

,  ,.n  2n+l 

(-1)  y 

La 

I2n+l| 

n“0 

M 

m 

v 

-  r-1  /2 

D_  ) 

La 

m  La  12 

rl 

n“0 

,  Nn  2n+2  2m-2n  .... 

(-1)  y  x  (93) 


The  stresses  and  displacements  can  thus  be  found  by  combining  equations 
(77)  through  (81)  and  (93): 


m 

°x 


6K,xy  +  2K  + 
1  8 

M 


_  r— »  2nH-l  , 

L  c"  I  Lil  <_1)'  <2,,fl>  2ny 


2n-l  2m-2n+l 


nr*2  n“0 
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f 


•i-%  a 


(-1)"  (2n+2)(2n+l)  y2n  x2m'2n 


(94) 


6K„xy  +  2KU  + 
M 


m 

Z  Cn>  Ym  (2I+1J  ('1>n  <2m-2n+l)(2m-2n)  y2n+1  x2^n-l  + 


nr*  2  n-0 

M  m 


Zr—  /2m+-n 

D"  L  2n+U  (-l)  (2m-2n) (2m-2n  -l )  y21*2  „2«-2"-2 


,  "  L-.  12n+l 

rrr- 1  n-0 


m 

rxy 


3Kly8  -  K3  -  3Kax» 
/2nrf  1' 


M  m 


Z  —  r-i  2nrt-l 

Cm  L  2n+l  (’1)n  <2nfl>(2m-2^1)  y2n  x2m'2n  - 


m-2  n-0 

H  m 


tfn+l, 
/2nrf  1 


Z-  r— ■  2rrrfl 

Dm  2^  2n+l  (_1)  <2rrf2)(2®-2n)  y2n+1  x2m-2n-l 


m-1  n-0 


E  u 
m  x 

l+vm 


[(1_vm)  (3xay-2y3)]  -  vmy3 

+  f.[-V  -  2|l-vm)y] 

r.  1  .1  -  r  .a 

['-Vmjy3]  +  Kn(-2vmx 

M 


(95) 


(96) 


kk+i  c„  1  r ]<-u". 

m-2  n*  J  2rrfl/ 

» 

I, ...  \  2n(2rri-l)  2n-l  2m-2rrf2  ,  . ,  ,  „  1 

ll  Vn>]  2m-2n+2  y  x  "  vm(2ra-2rrH)  y2n+1  x2m‘2nl  + 


(Continued) 


I  • 
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+ 1  5. 


m 


n“0 


(2n+2)(2n+l)  . 
2m-2n+l 


2n  2m-2n+l 

y  x 


vm(2m-2n) 


2n+2 

y 


2m-2n- 

x 


1 

■ 


(97) 


m 


«mU 


1+v. 


‘.[•V  -  I1-".)*’]  +  lt>(-'W‘)  +  Ks(-2v)  + 

*,  [(1-'v)0*y»  -  2*3)  -  V3]  +  Kn[2(i-v)y]  + 


M  m 

i  p  —  r'  1 

[2nrH\ 

['  1-Vn  x 

1+  Y.  c»  Z 

nr  2  n“0 

[2n+l/ 

(-1) 


\  (2m-2trH)(2m-2n)  2n+2  2m-2n-l 

2n+2  V  X 


y  x 

M 


/o  \  2n  2m-2n+l 
vm(2n+1)  y  x 


rrr  1  n“0 


2nrfl\ 

2n+l 


(-1)‘ 


N 


(2m-2n) (2m-2n-l)  ^2n+3  x2m-2n-2 


2n+3 


..  / rt  2n+l  2m-2n  . . . . 

vm(2n+2)  y  x  (98) 


A  basic  representative  element  of  fiber-reinforced  composites  under 
lateral  loading  is  depicted  in  Figure  34.  For  the  convenience  of  numerical 
calculation,  the  boundary  conditions  for  each  assigned  point  of  the  element 
are  listed  in  Table  I.  The  Cartesian  coordinates  of  all  points  are  listed 
in  Table  II.  The  normal  and  tangential  stresses  at  the  interface  are 
shown  in  Figure  35.  These  stresses  can  be  expressed  as 


c 


T 


n 


t 


•cp  +  cy  sir^cp  2Txy  sin  cp  cos  cp 

(99) 

sin  cp  cos  cp  +  TXy  cos^p-sin8cp 

(100) 
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*'*'  . --1  "IlililinTfl  i  MW 


f 


I 


1 


I 


> 


^mt 


Ox  coa^cp  +  Oy  sinacp  +  2<jXy  sin  cp  cos  cp 
|ay-oxJ  sin  cp  cos  cp  +  aXy|cos2qp-sinacp| 

u  cos  cp  +  v  sin  qp 
-u  sin  cp  +  v  cos  cp 


(101) 


At  poii.t  , 


i 

yi 


c  -  a  cos 


<Pi 


2  ~  a  ®ln  Cpi 


(102) 


Equations  which  satisfy  the  ten  types  of  boundary  conditions  are 
listed  in  this  section.  In  these  equations,  >01^ » •  •  •  ><*37  a>  •  •  •  * 

arl9J  identifies  the  19  unknowns  ,W0j  of  equa¬ 

tions  (88)  through  (92)  and  (94)  through  (98). 
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TYPE  1  EQUATIONS 


f  b  f] 
cx  c  "  fl  j  2 ^  | 


f  b  fl 
u  c-a,j,aij 


axy(c'a»2  ,ai) 

ml  b  ml 

°’xylc‘a*2  ,ail 

ml  b  ml 

-  u  |c"a*J»ui) 

v(c-a'l"i) 

v(ca,|«") 


0 


0 


0 


0 


0 


0 


(103) 

(104) 

(105) 

(106) 

(107) 

(108) 


TYPE  2  EQUATIONS 


m 


O 


y 


8  sin(  3o| 
ni]  k  . 

30/  ’  2 


+ 


0  (109) 


where 

j  “  1,...,1A 
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* 


A 


0  (110) 


0  (111) 
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-  |  Uy |c-a  cosjg]  ,  |  -  a  ainj^)  ,  or[j  - 

Uy[C‘a  COS(?o)  *  2  "  3  S±n(3o)  *  *i]  j  C°S(3o) 


(112) 


where 

j  =  1 , . . . , 14 


TYPE  3  EQUATIONS 


fib  f\ 

ml 

b 

ml 

<Mc>2  '  a*Q'l) 

C>2  - 

a*°'i) 

f  1 

b 

fl 

axyl 

c,2  " 

m  | 

b 

ml 

\ 

j 

^xyl 

°  ’2  ~ 

a'Q'i) 

1 

i 

f 

b 

f\ 

u 

X 

(c  *2  * 

a’«i  1 

\ 

i 

m 

b 

ml 

u 

X 

H  " 

a’ai! 

fib  t\ 

m 

b 

ml 

Uy|C>2  ’  a’«il 

-  u  | 

C’2  ■ 

a>ai] 

TYPE  4  EQUATIONS 


f  ..a  b  f 
CTxy  C'a+Jl0’2’ai  =  ° 


where 

j  -  i,  •  •  • ,9 


(113) 


(114) 


(115) 


(116) 


(117) 


(118) 


(119) 
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t 


where 


J  -  1, ... ,9 


TYPE  5  EQUATIONS 


f  ,  .  a  b  f 
u  c-a+tf-—  »or. 
y  *  2  i 


(120) 


i  -  1,2,3, 4 


f  b  f 
Txy  C‘2>0fi 


(121) 


f  b  f 


U  C’2'tti 


(122) 


ff  b  f 


Uy  C’2’"i 


(123) 


TYPE  6  EQUATIONS 


where 


J  -  1 . 9 


where 


j  -  1 . 9 


fib  ,  a  f 
CTxy  C’2_  *10’®!  "  0 


f  I  b  f 


u  rr  Jio’ai 


(124) 


(125) 


TYPE  7  EQUATIONS 


(126) 


k 


(127) 


where 

j  =  1,2,J,4 


TYPE  8  EQUATIONS 


TYPE  9  EQUATIONS 


where 

j  -  0,1... .,11 


where 

j  “  0,1 . 11 

TYPE  10  EQUATIONS 


where 

j  «  0,1,2 


m 
ot . 

i 


0 


m 


u 


k 


0 


m  c  _b  m 
CTxy 


0 


m 


u 

y 


12’  2’“i 


0 


m  ,e-a  _b  m 
axy[rT’2’Q'i 


0 


(128) 


(129) 


(180) 


(131 ) 


(132) 


(133) 
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The  total  number  of  equations  is  148,  but  the  total  number  of  unknowns 
to  be  determined  ,  af,  tfg  ,  •  •  •  .Qffg  )  is  19.  Therefore,  we 

can  use  the  method  of  Least  Squares  to  obtain  the  results.  When  solving 
the  148  simultaneous  equations  (103)  through  (134),  we  assumed  that  k  = 

-  1  . 

In  theory,  the  Point  Matching  method  was  used  to  solve  the  problem, 
and  the  computer  program  is  presented  in  Appendix  II.  The  finite  element 
method  was  eventually  used,  however,  since  it  was  found  that  the  Point- 
Matching  method  does  not  give  good  convergence  in  the  solution. 


APPENDIX  II 


S JOB »91 503-002.  TED  NEEF. 

*TAPE. SCR *02. OLD-OO. NEW-OO. 

RUNIS..... .240000) 

SET ( 0 ) 

LGO. 

PROGRAM  THREEDI INPUT  .OUTPUT ,  TAPE 5- INPUT • TAPE6-0UTPUT .TAPE  1 .TAPE2  I 
COMMON/1/  MSKI1000).  KS22  ( 80 ,240 . 3  )  .  I  VI  55 .4 .3  )  .XI  <  A ,  3 )  .X2 1 4 ,3 )  . 

1  X3I4.3),  K21I82)  .KMX  I  12 . 12  )  »CDA  (  72 .3 .55  >  . 

2  IN.  MM,  IP,  I  PL .  A21 1 82  >  .  XI32.2)  .  ICOI*. 3. 55) 

3  .  CENTO. 3. 55).  DZ 
DIMENSION  AS22I 30.2401 
EQUIVALENCE  IS( 19201  1.AS22) 

COMMON  /INPUT/  XNUI 2 ) .  £(2>  .XLAMI2)  •  G(2>  .  CMXI6.6.2)  .  DO 

DIMENSION  CI6.6.2)  ,  TEMPI  10) 

DIMENSION  Si  80 , 80, 5 ) 

EQUIVALENCE  I KS22 ( 19201 ) »S )  ,  (CMX.C) 

c 

C  THIS  PROGRAM.  WITH  ITS  ASSOCIATED  SUBROUTINES.  SOLVES  THE  THREE- 
C  DIMENSIONAL  STRESS  PROBLEM  BY  THE  METHOD  OF  FINITE  ELEMENTS. 

C 

C  MM  ••  HIGHEST  NODAL  POINT  INDEX  OF  TRIANGLES  IN  GROUND  PLANE 
C  IN  -  HIGHEST  INDEX  NUMBER  OF  BASE  TRIANGLE 

C  IP  -  HIGHEST  INDEX  OF  BASE  PLANE  NORMAL  TO  THE  Z  DIRECTION  HAVING 
C  THE  DISTANCE  IP*DZ  FROM  PLANE  Z  ■  0. 

MM  ■  32 
IN  *  55 
IP  ■  9 

READ  lOOO.CTEMPi I ) , I  *  1 , 1 0 ) 

1000  FORMAT  1 10A8) 

PRINT  1010.  (TEMPi I ) ,1-1,10) 

1010  FORMAT! 1H1.10A8) 

READ  1001.  DZ.Eil).  E(2>,  XNU(l).  XNU<2),  DO 

1001  F0RMATI6I5X.  E8.AJ) 

DO  1  I-l.MM 

1  READ  1005.  K .  XIK.1)  .  XIK.2) 

1005  FORMAT (I3.E12.4.E12.4) 

DO  10  1-1,2 

XLAM(l)  -  XNUI I ) *E  I  I  )  /  I ( 1 .♦XNUI  I ) )* 1 1 ,-2.«XNU( I ) )  ) 

G< I )  «  Ell)  /  12. *1 l.+XNUI  I  ) ) ) 

10  CONTINUE 

c 

C  DERIVE  C  MATRIX  CMX 
C 

DO  20  1-1 .6 
DO  20  J-1,6 
DO  20  K-1.2 
CMXII.J.K)  -  0. 

20  CONTINUE 
DO  25  1-1.2 

*  Cll.l.I)  -  XLAMII)+2.*GII) 

CI2.2.I)  -  Cll.l.I) 

CI3.3.I)  -  Cll.l.I) 

CIA. 4, I)  -  GUI 
CIS, 5, I)  -  CIA, A, I) 

CIS, 6, 1)  -  CIA. A, I) 

CI2.1.I)  -  XLAMII) 

CIS. 1,1)  -  CI2.1.I) 

CI3.2.I)  -  CI2.1.I) 

Cll.2,1 )  -  CI2.1.I) 

CI1.3.I)  -  CI2.1.I) 
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o  no  noo  oo 


C(2.3.I)  «  CI2.1.I) 

25  CONTINUE 

PRINT  OUT  THE  INPUT  QUANTITIES 

PRINT  1090.  ( XLAM( J ) . J« 1  .2  )  . (G( J) »J“1 .2 ) .( XNU t  J  >  »  J*1  *2  )  » 

1  (E(J).J-1.2I 

1090  FORMAT ( 1H1  .  5X .9HLAMBDA <  1 )  . 5X .9HLAMBDA< 2  > . 10X . 4HG< 1 1 . 10X.4HG 12 » t 

1  9X  »5HNU( 1 ) *9X  *5HNU  (  2  )  •  1 OX . AHE ( 1 ) « 10X  »4HE ( 2  )  /IX ,«< 3XE1 1 •*)/// 1 

PRINT  1091.  02.  DO 

1091  FORMAT  ( 1H  <  12X  »2HDZ  .  12X  .  2HDO/X  »2(3X*E11.4)///) 

DO  30  J  -  1.8 

11  -  4*0-1 )  +1 

12  -  11+3 

1092  FORMAT  (  1H  »4(  7X.  2HX (.I2.3H.1),  7X .2HX ( . 12 .3H.2  )  ) /X8 1 3XE1 1 .4  I / » 

3o  print  io92.(((ii.in.ii-n.i2).(ixm.i).xm*2)).n-ii»i2n 

PR 'NT  OUT  C  MATRICES.  1  PAGE 


PRINT  2001 

2000  FORMAT (2H  C.Il/I 
00  5001  I«l,2 
PRINT  2000,1 

00  5000  J-1,6 

5000  PRINT  2002.ICU.K.I )  »K«1.6> 

5001  PRINT  2003 

2001  FORMAT (1H1) 

2003  FORMAT!///) 

2002  FORMAT ( 1H  .6<E13#6.4X)> 

CALL  NOWAK 


NOW  HAVE  COEFFICIENTS  OF  UNKNOWN  DISPLACEMENT  MATRIX 
AS22 ( 80. 160 J  IS  FIRST  PARTITION  ROW 
KS22 (80X240 . 1 )  IS  USEO  FOR  ROWS  TWO  THROUGH  S 
KS22(80X240,2)  IS  GENERATED  IN  CHL30 


CALL  CHL30 

NOW  HAVE  UNKNOWN  DISPLACEMENTS  IN  SPIB2.9)  -  $11,1.11 


CALL  SIGMAS 

STOP 

END 
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SUBROUTINE  GETKMX(KJ.I) 

COMMON/1/  MSK(IOOO).  KS22  ( 80 . 240 . 3 )  *  I  V ( 55  .  4  .3 > .X 1 ( 4 , 3  )  ,  X  2  (  4  ,  3  )  , 

1  X  3  (  4  ♦  3  ) •  <21(821  »KMX  (  12.12), CDA(72»3»55)  , 

2  IN.  MM,  IP,  I  PL »  A21  ( 8?  )  ,  X(32.?)  ,  !C0(4,3,55) 

3  ,  CENTO. 3, 55)  ,  DZ 
DIMENSION  5(80.80.5) 

DIMENSION  AS22 ( 80 , 240 ) 

EQUIVALENCE  (S(  19201  1.AS22) 

EQUIVALENCE  ( KS22 ( 19201  )  »S  ) 

THIS  SUBROUTINE  DERIVES  THE  (12*  i  2)  K-^  '  MCES  FOR  USE 
IN  DERIVING  ROWS  OF  THE  BIG  VATRIX  <?? 

DIMENSION  DX  T ( 72  !  .  DX( 72) 

DIMENSION  JMX  I  1?  .  1  )  ,  AM*  !  2  > 

DIMENSION  XI  (  A  )  -  L  i  A  ( 4  )  ,  ZFi-  «l.  ASIRI4.4) 

COMMON  /INPUT/  XNU  ( 2  )  .  I  > ,  .XLAMOI  ,  O  ( 2 )  .  CMXI6.6.2)  ,  j 

REAL  KMX 

DATA  ( DXT ( J )  ,  J  =  1  »72)/0.»l.»16*0.»l. *16*0.  .1..0..0. ,1.  • 0 * , 0 • .1. ♦ 
1  9*0..1.»5*0.,1.,9*0.»1..0.,0.»1.»0./ 


DATA  (DX  (J)  »J=1.7c)/b*J«»l.»6*0.»l.»6*0.»l.»10*0.»l.»3*0.»l.» 
1  9*0.. 1. *10*0. , 1., 6*0. .1. ,2*0. .1* .3*0./ 


IF< I.LE.10) 

GO 

TO 

15 

IF( I.LE.12) 

GO 

TO 

10 

IF( I.LE.21  ) 

GO 

TO 

15 

IF( I.LE.24) 

GO 

TO 

10 

IF< I.LE.32) 

GO 

TO 

15 

10  IREG  -  2 
GO  TO  20 
15  IREG  -  1 
20  CONTINUE 


IREG  -  1.  FIBER 
IREG  ■  2.  RESIN 


DO  30  IBC*  2.4 

XI ( IBC  )  *  XI ( IBC  .KJ 1  -  Xl(’.KJ) 
ETA(IBC)  ■  X2IIBC.KJ)  -  X2I1.KJ) 
ZETA(IBC)  *  X3 ( IBC . K J )  -  L3I1.KJ) 
30  CONTINUE 


GET  INVERSE  OF  AMX 


ASTR(l.l)  ■ 

1 

2 

ASTR (2.1)  * 

1 

2 

ASTR (3,1)  « 

1 

2 

ASTR (4.1)  ■ 

1 

2 

ASTR ( 1.2)  ■ 
ASTR <  2.2)  * 
ASTRO. 2)  * 
ASTRO, 2)  « 
ASTRU.3)  « 


XI  (2)*(ETA(3)*ZETA(4)  -  ETA(4)*ZETA(3) ) 
-XI  ( 3)*(ETA(2)*2ETA(4)—  ETA (4  )  *ZETA( 2 ) ) 
♦  XI ( 4)*(ETA(2)*ZETA(3 )-  E TA ( 3 ) *ZETA ( 2 ) ) 
-(ETA(3)*ZETA(4)  -  ETA (4)*ZETA(3)  ) 
♦(ETA(2)*ZETA(4)  -  E TA ( 4 ) *ZE TA ( 2 )  ) 

- ( ET  A ( 2 ) *2ETA  ( 3 )  -  ET A ( 3 ) *ZE TA ( 2 )  ) 

(XI  <’ i*ZETA(4)  -  XI(4)*ZETA(3)) 

-IXI i2)*ZETA(4>  -  Xl(4)*ZETA(2)  ) 

+  ( X  I  (2  )*ZET  A( 3 )  -  X  I ( 3 ) *ZETA ( 2  ) ) 

-(XI ( 3 ) *ETA ( 4 )  -  XI ( 4 ) *E T A ( 3  )  ) 
♦(XI(2)*ETA(4)  -  XI (4)*ETA(2)  ) 

-(XI  (2)*ETA(3)  -  XI <3 )*ETA(2  )  ) 

0. 

I ETA( 3 ) *ZETA(4 )  -  ET A ( 4 ) *ZETA ( 3 )  ) 

-(XI (3>*ZETA<4)  -  XI (4)*ZETA(3) ) 

(XI ( 3 ) *ETA ( 4 )  -  XI (4>*ETA(3)  ) 

0. 
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ASTRI2.3)  *  -  (  ETA  <  2 )  *ZETA(4  )  -  ETA(4)*ZETA(2  >) 

ASTRO. 3)  *  (Xl  (2I*ZETA(4>  -  X  I ( 4 >*ZETA( 2 ) ) 

ASTR  (4.3)  ■  -(XI (2)*ETA(4)-XI(4)*ETA(2)  > 

ASTR ( 1.4)  =  0. 

ASTR (2.4)  =  ( ETA ( 2 ) *ZETA ( 3 )  -  ETA(3)»ZETA(2> > 

ASTR (3.4)  *  -(XI (2)*ZETA(3>  -  X  I (3 )*ZETA(2) ) 

ASTR (4.4)  »  (XI  (2)*ETA(3)  -  XI(3»*ETA(2)  ) 

C 

DELTA  »  X  I ( 2  )  *  (ETA ( 3 )*ZETA(4 I  -  ETA(4 )*ZETA( 3  ) ) 

1  +  X  I ( 3  )  *  (-ETA (2 )*ZETA( 4 )  +  ETA(4)«ZETA(2) ) 

2  +  X  I ( 4  )  *  ( ETA ( 2 ) *ZETA ( 3 )  -  ETA ( 3 ) «ZETA ( 2  )) 

DO  35  IQ  =  1  .4 

DO  35  IR  *  1.4 

35  ASTR ( 1 0. I R )  *  ASTR(IQ.IR)  /  DELTA 
DO  50  IQ  *1.4 

DO  50  IR  =1.4 

AMX(IQ.IR)  *  ASTR ( IQ. IR ) 

AMX ( I  Q+4, IR+4)  *  ASTR(IQ.IR) 

AMX ( IQ+8 . IR+8  )  *  ASTR ( 10 . 1 R ) 

AMX ( I Q+4 , I R  )  =  0. 

AMX ( IQ+8, IR  )  *  0. 

AMX( IQ+8.IR+4)  =  0. 

AMX (IQ  ,!R+4>  *  0. 

AMX ( IQ  , I R  +  8  )  «  0. 

AMX (IQ+4, IR+8 I  *  0. 

50  CONTINUE 
C 

C  NOW  AMX  CONTAINS  INVERSE  OF  A 
C 

V  *  ABS ( DELT  A  J  /  6. 

CALL  MXMULT (DX, AMX, KMX, 6, 12.12) 

CALL  MXMULT  (  CMX  (l.l.IREG). KMX, CDAd.KJ.  I »  .6.6,12) 

C 

C  NOW  HAVE  CDA(I) 

C  CDA  IS  (6  X  12).  SECOND  SUBSCRIPT  IS  TETRAHEDRON  NUMBER  KJ. 
C  THIRD  SUBSCRIPT  IS  BASE  TRIANGLE  INDEX. 

C 

DO  60  IQ  =1,4 
DO  60  IR  =1,4 
AMX(IQ.IR)  =ASTR ( IR, IQ) 

AMX ( I Q+4  » I R+4  )  = ASTR < I R , IQ ) 

60  AMX( IQ+8, IR+8)  =ASTR(IR,IQ> 

C 

C  AMX  NOW  CONTAINS  TRANSPOSE  OF  INVERSE  OF  A 
C 

CALL  MXMULT ( AMX, DXT.JMX.12  > 12,6  ) 

CALL  MXMULT ( JMX  .CDA ( 1 ,K J, I ) »KMX ,12,6.12) 

CALL  MXCON(KMX.KMX,V,12,12) 

C 

C  K  MATRIX  NOW  IN  KMX(J.K)  ,  J*1.12  »  K*1.12 

C 

40  RETURN 
END 
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SUBROUTINE  NOWAK. 

COMMON/1/  MSK(IOOO).  KS22 ( 80 ,240 , 3 )  .  I V( 55  ,4 ,3 > ,X 1  ( 4 , 3 ) ,X2 ( 4 ,3 ) » 

1  X 3 ( 4 , 3 ) »  K2 1 ( 82  )  ,KMX( 12 ,12) »CDA(72»3,55 )  , 

COMMON/1/  MSM1000).  KS22  (  80 ,240, 3  )  ,  I V  (  55  ,4  ,3  )  ,X  1  ( 4 , 3  )  ,X2  <  4 , 3 )  , 

1  X  3 ( 4 , 3  )  ,  K2 1 ( 82  )  *KMX (12,12) ,CDA( 72,3,55)  , 

2  IN.  MM,  IP.  I  PL ,  A21I82)  ,  X(32,2>  ,  ICO(4,3,55) 

3  ,  CENTO, 3. 55) .  OZ 

DIMENSION  AS22 ( 80,240 )  ,KX ( 90 ) ,PAZ ( 3  ) 

EQUIVALENCE  (S(  19231  ) .AS22  ) 

COMMON  /INPUT/  XNUI2),  E(2)  ,XLAM(2)  .  G<2)  ,  CMX(6,6,2)  ,  DO 
EQUIVALENCE  (KS22I19201)  ,QS22 ) 

DIMENSION  S ( 80 , 80,5 ) 

DIMENSION  QS22I 12800)  ,  PP(1000) 

EQUIVALENCE  (KS22( 19201 ) «S) 

DIMENSION  IM( 35.3)  .  LMI35.3) 

DIMENSION  1X0(3.21 
REAL  KS22.  K21  .  KMX 
C 

C  THIS  SUBROUTINE  DERIVES  ALL  NON-ZERO  COEFFICIENTS  FOR  THE  UNKNOWN 
C  DISPLACEMENT  MATRIX.  THE  DISPLACEMENT  MATRIX  IS  (720X720).  AND  IS 
C  DIVIDED  INTO  9  ROWS  OF  SUBMATRICES  WHICH  HAVE  (80X80)  ELEMENTS  EACH. 
C  AT  MOST,  THREE  OF  THESE  SUBMATRICES  CONTAIN  NON-ZERO  ELEMENTS.  THUS. 
C  DIVIDED  INTO  9  ROWS  OF  SUBMATRICES  WHICH  HAVE  (80X80)  ELEMENTS  EACH. 
C  AT  MOST.  THREE  OF  THESE  SUBMATRICES  CONTAIN  NON-ZERO  ELEMENTS.  THUS, 
C  ONLY  THESE  THREE  NON-ZERO  SUBMATRICES  ARE  DEVELOPED. 

C  SUBROUTINE  NOWAK  PUTS  THE  TWO  SUBMATRICES  OF  INTEREST  FOR  THE  FIRST 
C  PARTITION  ROW  INTO  AS22 ( 80X160 ) .  FOR  ROWS  TWO  THROUGH  EIGHT.  THE 
C  NON-ZERO  SUBMATRICES  REPEAT.  THESE  THREE  MATRICES  ARE  STORED  IN 
C  KS22 ( 80X240, 1 ) • 

C  THE  TWO  PARTITION  ELEMENTS  IN  ROW  NINE  ARE  GENERATED  IN 
C  SUBROUTINE  CHL3D  WHEN  NEEDED  ON  THE  NINTH  PASS. 

C 

C 

C  THE  INDEPENDENT  TERMS  CORRESPONDING  TO  ROW  ONE  WILL  BE  IN  A21. 

C  THE  INDEPENDENT  TERMS  FOR  ROWS  TWO  THRU  NINE  WILL  BE  IN  K2K80.1) 

C 

C 

c 

C  XI  1. 1)  -  X  COORDINATE  OF  ITH  NODAL  POINT 
C  X( I  ,2)  -  Y  COORDINATE  OF  ITH  NODAL  POINT 
C  DZ  -  DISTANCE  BETWEEN  TWO  CONSECUTIVE  BASE  PLANES 
C  E(l)  -  MODULUS  OF  ELASTICITY  OF  MEDIUM  1 (FIBER) 

C  E ( 2 )  -  MODULUS  OF  ELASTICITY  OF  MEDIUM  2 (RESIN) 

C  XNU(l)  -  POISSONS  RATIO  OF  FIBER 
C  XNU ( 2 )  -  POISSONS  RATIO  OF  RESIN 
C  I  CO  (  4 , 3  )  -  ZERO  OR  ONE 

C  IKO( 3,2  )  -  ZERO  OR  ONE 

C  DO  -  DISPLACEMENT  IN  X  DIRECTION  ATPLANE  X-  ( SORT ( 3. ) /2 . >*B 
C 

PYF  *  0. 

PZF  *  0. 

PXS  *  0. 

PZS  *  0. 

C  P  IS  FORCE. 

C  F  IS  FRONT  •  S  IS  SIDE,  X,Y,Z  ARE  DIRECTIONS 
C 

PAZ(l)  *  0. 

PAZ ( 2 )  ■  0. 

PAZ ( 3  )  -  0. 

C 

C  PAZ( 1  )  IS  X  COMPONENT  OF  FORCE  AT  PLANE  Z*0(B0UNDARY  LOAD) 

C  PAZ ( 2  )  IS  Y  COMPONENT  OF  FORCE  AT  PLANE  Z«0( BOUNDARY  LOAD) 
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PAZ ( 3 )  IS  Z  COMPONENT  OF  FORCE  AT  PLANE  Z-0(BOUNDARY  LOAD) 

FOR  PAZ ( 1 )  -  PAZ (2)  -  PAZ ( 3 ) *  0  WE  HAVE  FREE  SURFACE  CONDITIONS. 

IX  «  0 
JX  «  0 

DO  103  M  ■  1 .MM 
IFIM.LT. 8)  GO  TO  100 

IF(M0D(M»7 ) . EQ.O.OR.MOD(M-1 .7 )  .EO.O )  GO  TO  101 
DO  105  MNl-1.3 
IX  ■  IX  ♦  1 
JX  «  JX  +  1 

MSK(IX)  -  3«<M-1)  ♦  MN1 
PP(IX)  ■  PAZ(MNl) 

105  KX(JX)  -  IX 
GO  TO  103 

100  IF(M.E0.1.0R.M.EQ.7)  GO  TO  102 
IX  ■  IX  +  1 

MSK ( IX ) *  3«(M-1)  ♦  2 

PP(IX)  «  PYF 

IX  ■  IX  +  1 

MSK ( IX )  -  3* ( M-l )+3 

PP(rX)  «  PZF 

GO  TO  103 

102  IX  «  IX+1 

MSK (IX)  *  3*M 
PP(IX)  -  PZS 
GO  TO  103 

101  IX  ■  IX  ♦  1 

MSK (IX)*  3  *  (M-l)  ♦  1 
PP(IX)  »  PXS 
IX  »  IX  +  1 

MSK ( IX )  ■  3» (M-l )+3 
PP(IX)  »  PZS 

103  CONTINUE 
IXMAX  ■  IX 
JXMAX  »  JX 

DO  104  N* 1 » 1 1 
DO  109  IX  ■  It  IXMAX 
MSK  ( N*  I XMAX+ 1 X  )  »MSK  (IX)  +  N  *  3  •  MM 
109  PP( N* I XMAX+IX )  ■  PP(IX) 

DO  108  JX  «  1. JXMAX 
KXJX»KX( JX) 

108  PP(KXJX+N* IXMAX)  ■  0. 

104  CONTINUE 

DO  169  I  ■  1.  IN 
IFF  *  1 

IFIMOD ( 1.2) .EO.O )  IFF  ■  -1 
IF(MOD<I.12).EO.O)GO  TO  2 
I  El  2  »  0 
GO  TO  4 
2  I E 1 2  *  1 

4  IM(l.l)  «  1  +1/2  +  U-IE12J/12 

IKO(ltl)  «  0 
IZ  -  1/2 

IF(MOD( IZ.2I.EQ.0)  GO  TO  6 
IF(I.GT.48)  GO  TO  12 
I M ( 2 *  1 )  »  IM(I.I)  +IFF  +  7 
IK0I2.1)-  0 
GO  TO  7 

12  IM( 3.1 )  «  IM( 1.1 )+IFF 
IK013.1)  »  0 
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GO  TO  7 

6  IF< I.GT.48)  GO  TO  11 
IM(2.1)  ■  IM(1.1)  IFF 
I  KO ( 2  » 1 )  ■  0 

GO  TO  7 

11  I M ( 3*1  )  «  IM( 1 *1 )  ♦  IFF  +  7 
IKO ( 3  » 1 )  *  0 

7  I F (  I.GT.48)  GO  TO  13 
IM( 3 » 1  )  *  I M < 1 • 1 )  +  7 
IKO ( 3  *  1 )  *  0 

GO  TO  30 

13  IM( 2  » 1  )  -  IM(1.1)  7 

I KO ( 2  ♦  1 )  *  0 

30  1F(  IM(1,1).GT.32)  GO  TO  31 
GO  TO  32 

31  IM< 1.1 )  *  64  -  IMI1.1 > 

IKO (1.1)  «  1 

32  IF( IM(2.1 J.GT.32)  GO  TO  33 
GO  TO  34 

33  IM(2.1)  *  64  -  IM(2.1) 

IKO (2. II  *  1 

34  IF( IM(3.1).GT*32)  GO  TO  35 
GO  TO  36 

35  IM( 3 . 1  )  ■  64  -  IM( 3 . 1  ) 

IKO (3.1  I  «  1 

36  DO  20  K  *  1.3 

DETERMINE  ALPHA.  BETA,  GAMMA,  FOR  EACH  NODE  IN  EACH  BASE  TRIANGLE  I 
DETERMINE  ALPHA,  BETA,  GAMMA.  FOR  EACH  NODE  IN  EACH  BASE  TRIANGLE  I 


K* 1  CORRESPONDS  TO  POINT  ALPHA 
K-2  CORRESPONDS  TO  POINT  BETA 
K*3  CORt-TcSPCNDS  TO  POINT  GAMMA 

I  MM  =  MOD( IMIK.1 ) ,21 > 

IF ( IMM.EO*0)  I  MM* 21 

I F ( MOD ( I  MM , 7 ) . EO . 0 )  GO  TO  8 
IE7*  0 
GO  TO  10 
8  IE7  *  1 

10  IBAR  «  ( IMM-IE7)  /  7 
ID2  *  MOD ( IMM.2 ) 

IH1  =  MOD ( IBAR, 3 )  -1 
IF(IHl)  14.15,16 

14  I F ( ID2.EQ.0)  GO  TO  17 
I M ( 3.2)  *  IM(K.l) 

IKO ( 3 ,2  )  *  IKO(K.l) 

GO  TO  20 

17  IM( 1 » 2  )  *  IM(K.l) 

I  KO  ( 1  '<  2  )  -  IKO(K.l) 

GO  TO  20 

15  I F ( ID2.EO.O)  GOTO  18 
IM ( 3 » 2  )  =  IM (K , 1 ) 

I KO (3.2)  *  IKO(K.l) 

GO  TO  20 

18  I M ( 2,2)  *  IM(K.l) 

I KO ( 2 , 2  )  *  IKO(K.l) 

GO  TO  20 

16  IF  ( ID2.EQ.0)  GO  TO  19 
IM( 1,2)  »  IM(K.l) 
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c 

c 

c 

c 

c 
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IK0(1.2)  -  I KO ( K  » 1 ) 

GO  TO  20 

19  IM ( 2 * 2 )  -  I M ( K • 1 ) 

I KO ( 2  »  2 )  =  IKO(K.l) 

20  CONTINUE 

IFCMODI I.2I.E0.1 )  GO  TO  27 

UC  *  2 

MX  «  3 

MM1  *  MM 

MM2  «  0 

XM31  *  DZ 

XM32*  0. 

GO  TO  26 
27  LX  =  3 
MX  *  2 
MM1  *  0 
MM2  *  MM 
XM31  =  0. 

XM32»  DZ 

26  IV<  1*1*1)  =  I M ( 1  *  1 ) 

I V ( 1*2*1)  =  IM(LX.l) 

I  V(  I  »  3 » 1 )  *  IM(MX.l) 

I V( 1,4,1)  *  I M ( 2 ,2 )  +  MM 
I V( 1,1,2)  ■  IM( 1 .2 ) 

IV(  1,2.2)  >■  IMILX.2)  ♦  MM1 
IV<  1,3,2)  -  I M (MX , 2 )  +  MM2 

I  VC  I  ,4,2)  ■  IMI3.2)  +  MM 
I V ( 1*1*3)  *  IM(I.I)  +  MM 
I  VC  1,2,3)  =  I M ( MX , 1 )  +  MM 
I  VC  1,3.3)  =  IM(LX.I)  +  MM 
I V ( I ,4,3)  =  I M ( 1 ,2 ) 

MATRIX  IV  RELATES  THE  NODES  TO  THE  TETRAHEDRONS 

FIRST  INDEX  OF  IV  IS  THE  BASE  TRIANGLE  INDEX  I 

SECOND  INDEX  IS  NODE  POSITION  WITHIN  TETRAHEDRON  I 

THIRD  INDEX  IS  THE  NUMBER  OF  THE  TETRAHEDRON  ABOVE  BASE  TRIANGLE  I 

THIRD  INDEX  IS  THE  NUMBER  OF  THE  TETRAHEDRON  ABOVE  BASE  TRIANGLE  I 


X  3 ( 1  ,1  ) 

= 

0. 

« 

X3I2.1  ) 

= 

X  3  (  1,1  ) 

X3( 3, 1  ) 

= 

X3  (  1  .1  ) 

X  3  1  4 . 1  ) 

= 

X  3  C  1  , 1  )  +  DZ 

X  3  C 1,2) 

= 

X  3  (  1  . 1  ) 

X3I2.2) 

= 

XM3 1 

i 

X3( 3,2  ) 

= 

XM32 

X3 ( 4 , 2  ) 

= 

X3 (4,1 ) 

* 

X  3  C 1 ,3) 

= 

X3 ( 4 , 1 ) 

X3 ( 2  » 3 ) 

= 

X  3  (  4 , 1  ) 

% 

X3I3.3) 

= 

X3 ( 4 , 1 ) 

X  3  (  4  •  3  ) 

= 

X  3  C  1,1 ) 

i 

XIII.JI  IS 

THE  X  COORDINATE 

OF 

THF 

ITH 

NODE 

OK 

THE 

JTH 

TETRAHEDRON 

X?(I,J>  IS 

THE  Y  COORDINATE 

OF 

THE 

ITH 

NODE 

OF 

THE 

JTH 

TETRAHEDRON 

t  y 

X3(I.J)  IS 

THE  Z  COORDINATE 

OF 

THE 

ITH 

NODE 

OF 

THE 

JTH 

TETRAHEDRON 

ABOVE  THE  CURRENT  BASE  TRIANGLE 


ICO ( 1  , 1  .  I  ) 
I  CO  I  2 . 1 .  I  ) 
ICO! 3.1 . 1  ) 
ICO ( 4 . 1  . 1  ) 
I  CO  ( 1 , 2  .  I  ) 


IKO(l.l) 
IKOILX.l  ) 
IK0IMX.1  ) 
I KO I  2 , 2 ) 

I  KO (1,2) 
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I  CO (2*2.1)  *  IXOILX.2) 

ICOI3-2.I)  =  IXOIMX.2) 

1  CO (4.2.1 )  =  1X0(3. 2) 

I  CO ( 1.3.11  =  t KOI  1,1) 

ICO(2.3.I)  *  IXO(MX.l) 

icon. 3. n  *  ixoax.i) 

I  CO (4, 3.1)  =  1X0(1.21 
DO  50  IC1  *  1.3 
DO  50  IC2  «  1.3 
50  CENT ( IC1.IC2.I )  =  0. 

DO  98  XJ  =  1,3 

Xl(l.XJ)  =  X ( MOD ( IV( I.l.KJ) »MM)»1)*(-1)**IC0(1»XJ»I ) 

X  2  I  1  ,  X  J  )  =  X (MODI  I V(  r . 1 ,XJ)  ,MM)  ,2  )  * ( - 1 )** ICO( 1 .XJ » I ) 

X1I2.XJ)  =  X ( MOD ( 1 V (  I «2<KJ) »MM  ) ,1 )*(-l)**IC0(2»XJ» I  ) 

X2I2.XJ)  =  XIMODI I V (  I ,2.XJ) ,MM) .2 ) * ( - 1 ) **  I  CO  I  2 ,X J  ,  l ) 

X1I3.XJ)  =  XIMODI IVI I .3.XJ) ,MM) ,1 >*(-1 >**IC0(3.XJ,I ) 

X2I3.XJ)  *  X (MODI IVI I .3.XJ ) ,MM) ,2 ! *(-l )**ICO( 3 .XJ, I ) 

X1I4.KJ)  =  X (MOD  I  I  VI  I  ,4,XJ)  ,MM)  ,1  )* (-1 )** IC0I4.XJ,  I  ) 

X2I4.XJ)  =  XIMODI  IVI  I ,4,XJ)  ,MM) ,2)*("1 >**IC0(4.XJ, 1  ) 

C 

C  CENT  I  I  , J ,K )  IS  THE  ITH  CENTROID  COMPONENT  OF  THE  JTH  TETRAHEDRON, 
C  ABOVE  THE  XTH  TRIANGLE 

C 

C  CALCULATE  THE  CENTROIDS 

C 

IC1  =  XJ 
DO  55  I C 3 * 1  .4 

CENT  I  1  ,  IC1 , 1  )  »  CENT  I  1  «  I C 1  ,  I  )  +  XKIC3.IC1) 

CENT  I  2 ,  IC 1 , 1 )  *  CENTI2.IC1.I)  ♦  X2IIC3.IC1) 

55  CENT  I  3 , IC1 » I )  =  CENT  I  3 , I C 1  ,  I  )  *  X3IIC3.IC1) 

DO  56  IC4  =  1,3 

CENTIIC4.ICl.il  »  CENT! IC4.IC1.I)  /  4. 

56  CONTINUE 
C 

CALL  GETXMX I  X J , I ) 

421  DO  99  X  *  1.2 

DO  799  NZ  *  1.12 

ISA1  =  1 

Ml  *  M0DINZ.4) 

IF(Ml.EQ.O)  GO  TO  700 
IE4  =  0 
GO  TO  800 
700  Ml  =  4 
IE4  =  1 

800  M2  =  1  +  (NZ  -  IE4)  /  4 

IFl ICOIMl.XJ.I ).E0.1 )  GO  TO  799 
I  PL  =  3*1  IV  l  I  .Ml  .XJJ-MX-l  )*MM-1  >+M2 
I E I IPL.LE.ll I  GO  TO  300 
I  PL  =  I  PL -10* I  I  PL/ 96  + 1 ) 

300  CONTINUE 

IFl 3*1 IVI I.M1.XJJ+IX-1 )*MM-1 J+M2-MSXI IPL ) ) 301 , 302 ,303 
302  IFl IPL. GT. 160)  GO  TO  799 
III  =  IPL 
GO  TO  327 

301  IFIIPL.LT. 3)  GO  TO  799 
IMP  =  IPL  -  1 

306  I F I  3  *  I IVI I »M1,KJ)  +  (X-1 ) *MM-1 ) +M2-MSX I IMR ) )  304,305.799 
305  IFl IMR. GT. 160)  GO  TO  799 
III  *  IMR 
GO  TO  327 
304  IMR  =  IMR  -  1 
GO  TO  306 
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303  I MR  ■  IPL  ♦  1 

307  IF(3«<  IV(  I.Ml.KJl-MK-l  )*MM-1  )+M2-MSK<  IMR)  )  799.305.310 
310  IMR  ■  !MR*1 
GO  TO  307 

327  !F< 1 1 1 .LT.B1 )  GO  TO  391 
GO  TO  392 

391  A21  ( 1 1 1 )  «  A21lim+PP  (111) 

GO  TO  395 

392  K21UI1-  80)  ■  K2HII1-  80)  ♦  PP (III) 

395  DO  81  KA  «  1.3 
00  82  KB  «  1.4 

IPL  ■  3*( I V ( I .KB.KJ)  +  <K-i )*MM-1 >+KA 
I F ( IPL.LE.11)  GO  TO  396 
IPL  ■  IPL-10*< IPL/96+1) 

396  CONTINUE 

IFf  3*  (  I V (  I,K8,KJ)+(K-1  )*MM- 1 > +KA-MSK <  IPL)  )  349.350.351 

350  II 2  *  IPL 
ISA2  ■  ISA1 

IFOCA.LT. 3)  1 5A2  •  1SA1* ( 1-2* I  CO  I KB.KJ » I > ) 

OX  4X  487 

IF(KA.LT.3)  I SA2  *  ISA  1*11-2*1 CO (KBtKJ.I ) ) 

GO  TO  376 

349  IF(IPL.EQ.l)  GO  TO  374 
IF ( IPL.EQ.21  GO  TO  82 
IMR  -  IPL  -  1 

353  IF(3*( I V ( I.K8.KJ)+(K-1 ) *MM-1 ) +KA-MSK ( IMR ) )  354,352,375 

354  IMR  «  IMR  -  1 
GO  TO  353 

352  1 1 2  •  IMR 
I SA2  «  ISA1 

IFfKA.LT. 31  ISA2  ■  I SA1*( 1 -2* ICO ( KB .KJ ,  I  ) ) 

GO  TO  376 

351  IMR  -  IPL  +  1 

358  IF<3«(IV(I,KB,KJ)-MK-l)*MM-l)«-|f*-MSKUMR))  375.352,364 
364  IMR  -  IMR  ♦  1 
GO  TO  358 

376  IFf  II1.LT.81)  GO  TO  517 
GO  TO  519 

517  AS22 1  1 1 1  «  1 12  )  *  AS22  (  1 1 1 ,  I  1 2  )  +  KMX  (NZ.4*  (  KA-1  )  -t-KB )  *  ISA2 
GO  TO  82 

519  KS22I III-  80.112  .1)  ■  KS22III1-  80,112  ,1)  ♦ 

1  KMX(NZ»4*IKA-1 )*KB)*ISA2 

GO  TO  82 

375  MLD  «  M0DI3*( IVf I ,KB»KJ)-1 )+KA ,3*MM) 

IF  f  MLO.GT . 19 .OR. MLD.EQ.2 )  GO  TO  82 
374  IFf  II1.LT.81 )  GO  TO  518 
GO  TO  520 

518  A2KII1)  ■  A21 1  III ) -KMX  (NZ.4*  (KA-1  )*KB)*DO 
GO  TO  82 

520  K2KII1-80)  «  K2KII1-80)  -  KMX  (NZ  .4*  (KA-1  )+KB  )  «D0 
82  CONTINUE 

81  CONTINUE 
799  CONTINUE 
99  CONTINUE 
98  CONTINUE 
169  CONTINUE 
RETURN 
END 
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SUBROUTINE  CHL30 
C 

COMMON/1/  MSKdOOO).  <S22 ( 80 . 240. 3 ) . I  V ( 55 .4 . 3 ) .X 1 ( 4 , 3 )  .  X2 ( 4 . 3 ) • 

1  X3(4,3>,  <21(82)  .KMX ( 12 , 12 > .CDA I  72 .3 .55 )  . 

2  IN.  MM.  IP.  I  PL .  A21I82)  »  XI32.2)  .  ICOI4.3.55) 

1  X3I4.3).  <21(821  .KMXI 12.12) .CDAI72. 3.55 )  . 

2  IN.  MM.  IP.  I  PL .  A21I82)  •  X ( 3 2.2)  »  ICOI4.3.55) 

3  .  CENT (3.3.55).  DZ 
DIMENSION  JRCI32) 

DIMENSION  <9(80.160) 

EQUIVALENCE  <S.<9) 

DIMENSION  AS22( 80.240) 

EQUIVALENCE  (S(  19201  ) . AS22  ) 

DIMENSION  DELXI80.9)  •  GSAV(80,9) 

EQUIVALENCE  ( AS22  .2) 

DIMENSION  2(12800) 

DIMENSION  S( 80.80.5 ) 

DIMENSION  AS ( 80 . 80 ) 

EQUIVALENCE  (<S22( 19201 ) .S) 

C  THIS  ROUTINE  SOLVES  SU«G  .  WHERE  S  IS  A  TRI-DIAGONAL  MATRIX  IN 
C  SUBMATRICES.  WITH  ELEMENTS  OF  ORDER  N. 

C  S  IS  <NOWN 

C  SP  IS  A  VECTOR  OF  DIMENSION  (NXM I  WHERE  M  IS  THE  NUMBER  OF  DIVISIONS  OF  S 

C  CIS  WRITTEN  ONTO  TAPE  AFTER  DERIVATION  ON  THE  FORWARD  PASS. 

C  SP  IS  A  VECTOR  OF  DIMENSION  (NXM)  WHERE  M  IS  THE  NUMBER  OF  DIVISIONS  OF  S 

V.  CIS  WRITTEN  ONTO  TAPE  AFTER  DERIVATION  ON  THE  FORWARD  PASS. 

C  AND  READ  BAC<  IN  ON  THE  BAC<SWEEP 
C  S(lrU’)  INITIALLY  CONTAINS  S(I.I-l) 

C  Stl.lt?)  INITIALLY  CONTAINS  S(I.I  ) 

C  S ( i  .1.4*  INITIALLY  CONTAINS  SII.I  +  1) 

C  SP  CORRFSPONDS  TO  P  IN  THE  WRITEUP  BY  GATEWOOD  ON  THE  FORWARD  PASS. 

C  ON  IHT  P4CKSWEEP.  IT  CORRESPONDS  TO  U. 

Li  SION  SP ( 80 .  9) 

0Il  SION  C ( 6400 ) 

EOL  .  .ENCE  ( S ( 25601 ) ,C  ) 

EQU*.  -ENCE  (  AS22I12801 ) .SP) 

DIMENSION  Q( 12800) 

EQUIVALENCE  (S.Q) 

REAL  KS22 .<2 1 
REAL  <9 

DIMENSION  G ( 82  ) 

DIMENSION  ND ( 80 ) 

DIMENSION  SQI80.9) 

DATA  ( NO ( I ) . I *1 .80)  /  3.  2.  3.  2.  3.  2.  3.  2.  3.  2.  3.  3.  1.  3. 

1  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  3.  1.  3.  1. 

2  2.  3.  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  3.  1.  3.  1.  2.  3, 

3  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  2.  3.  1.  3.  1.  3.  1.  2.  3.  1. 

4  2.  3.  1.  2.  3  / 

DATA  I JRCd ) .  1-1  .32)  /  1.  3.  5.  7.  9.  11  ,  12.  14,  17,  20.  23. 

1  26,  29.  31.  33.  36,  39.  42,  45,  48,  50,  52,  55.  58.  61,  64, 

2  67,  69.  71,  74,  77,  80  / 

XLIMIT  *  l.E-8 

REWIND  1 
REWIND  2 
N  *  80 
N2  *  160 
<REM  «  80 
<REM2  *  16U 
M  .  9 

DO  30  ICYCLE  *1.M 
I F ( ICYCLE. GT . 1 )  GO  TO  12 
WRITE! 2)  <Z(  I  ), 1-1,12800) 
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DO  11  1-1.80 
G< I)  »  A21< I) 

GSAVI l.ICYCLE)  -  G<  It 

DO  11  J-1.80 

S(  I . J.2 )  -  AS22 ( I .J) 

11  S( I ,J.4)  «  AS22 ( 1  •  J+80 ) 

GO  TO  17 

12  I F(  I  CYCLE  .EO.M )  GO  TO  14 
DO  13  1-1.80 

G(  I)  -  K.  2 1  (  I  ) 

GSAVI l.ICYCLE)  »  G(I) 

DO  13  J-1.80 

S( I.J.l)  -  KS22I I • J . 1 ) 

S( I  .J.2)  -  KS22  t I  .J+80,1 ) 

13  S I  I . J .4 )  -  KS22I I . J+160,1 > 

GO  TO  17 

14  INEW  -  0 
JJRC  -  1 

DO  480  J-1.80 
GUI  -  K21 1 J  ) 

GSAVI J.ICYCLE)  -  G ( J ) 

IF  ( JRCUJRO.LT. J)  JJRC  -  JJRC  +  1 
I F ( J-JRC< JJRC) )  460.465.460 
460  INEW  -  INEW  +  1 

MSK  I  INEW+720  )  «  MSKIJ+720) 

MSK  (  I NEW+800  )  *  MSKIJ+800) 

465  DO  480  1-1.80 

SI  I .J.l)  «  KS22I I.J.l) 

I F ( J-JRC I JJRC ) )  470.475,470 
470  SII.J.2)  «  KS22I  I.J+80.1)  +  ICS22  (  I  .J+160.1  ) 

GO  TO  480 

475  SII.J.2)  -  KS22I I .J+80.1 ) 

480  CONTINUE 

WRITE  (2)  (01  I ).  1-1 .12800) 

DO  490  J-1,48 
MSKIJ+768)  -  MSKIJ+800) 

490  CONTINUE 
17  CONTINUE 
1  K1  «  N 
K2  -  N 
K  3  «  N2 
K4  «  N 
4  CONTINUE 

IF! ICYCLE.EO. 1 )  GO  TO  10 

CALL  MXMULT (Sll.1,1)  ,  511,1,5)  .  511,1,3)  .K1.IC2.K1) 

SIl.1,5)  CONTAINS  C  FROM  LAST  CYCLE 

CALL  MXSUBISI 1.1.2)  .  511,1,3)  ,  SIl.1,2)  .Kl.Kl) 

10  CONTINUE 

B(I.I)  NOW  IN  SIl.1,2) 

CALL  I NVERT (SIl.1,2)  .  K1  .  X3.  XLIMIT  .  FLAG) 

IF  (FLAG.NE.O.)  GO  TO  500 

INVERSE  OF  B I  I , I )  NOW  IN  511,1,2) 

IF  ( ICYCLE.EO. 1 )  GO  TO  20 

CALL  MXMULT  (511,1,1)  .SP I  1  .  I C  YCLE- 1 )  ,  SU.1.3)  ,Kl,  K2  ,  1) 
CALL  MXSUB  (G  .  SU.1.3).  G.Kl,  1) 

20  CONTINUE 
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CALL  MXMULT(S( 1.1.2)  .  G  .SP(1.  ICVCLE  )  .  K'  .  K1  .  1) 

if  < icycle.ge.mj  r.o  to  35 

CALL  MXMULT  (5(1.1. 2)  •  5(1. 1,4)  .  Sll.1.5).  XI  .  XI  .X4) 
NSO  -  N*X4 

WRITE  (1)  (C(  J) .J-l.NSO) 

30  CONTINUE 

35  CONTINUE 

NOW  IN  BACXSWEEP.  SOLVING  FOR  U 

DO  60  I  «  2.M 

JCYCLE  «  M-I+l 

IF( JCYCLE.LT. M-1J  GO  TO  36 

XI  -  XREM 

GO  TO  37 

36  XI  •  N 

37  CONTINUE 
NSO  *  X 1*N 

IF ( JCYCLE. EQ.M-1  )  GO  TO  41 
BACXSPACE  1 
41  CONTINUE 
BACXSPACE  1 

READ  (1)  (C(J) .J-l.NSO) 

U(M)«SP(M)  ,  CONSIDER  FIRST  (fc-l)TH  CYCLE 

CALL  MXMULT (S( 1.1.  5  I  .SPI 1 . JCYCLE+1 ) »S( 1 .  1  .  1  )  .N.Xl ,  1) 
CALL  MXSUB(SP( 1. JCYCLE).  S(l.l.l)  ,SP ( 1 .JCYCLE ) .  N  ,  1) 

60  CONTINUE 

U(NS.l)  NOW  STORED  IN  SP(N.I)  .  I-l.H 

REWIND  2 
DO  400  IC  -  1.9 
IFUC-2)  300,320,320 
300  READ  (2)  ( 0( I  ) . I -1 . 12800 ) 

CALL  MXMULT ( S. SP ( 1. 1  I. DELX(  1.1  ) .80, 160.1) 

GO  TO  400 

320  I F( IC.EQ.M)  GO  TO  350 
310  DO  315  1-1.80 
DO  315  J-1.240 

315  XS22U.J.2)  «  XS22  (  I  •  J,  1  ) 

CALL  MXMULT (S.SP(1«IC-1 ) .DELX ( 1 , 1 C ) .80 » 240 . 1 ) 

GO  TO  400 

350  READ  (2)  (Q( I ) .1-1,12800) 

CALL  MXMULT  (0,SP(1. 8  )  ,DELX(  1,10,80,160,1  ) 

400  CONTINUE 
PRINT  '001 

4001  FORMAT (iH1.5X,3HROW,l IX, 4HDELX.12X.3HX21  ) 

00  450  1-1,9 
DO  450  J-1.80 
X  -  80* ( I -1 ) ♦ J 

I F ( MOD (X.50J.EQ.0)  PRINT  4001 
450  PRINT  4000.X»DELX(J»,).GSAV(J»I ) 

4000  FORMAT(6X»I3.2(2X.E13,6)  ) 

RETURN 

500  CONTINUE 

PRINT  1000  ,  ICYCLE 

1000  FORMAT  ( 31H1C0ULD  NOT  INVERT  MATRIX  IN  ROW. 12) 

STOP 

END 
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SUBROUTINE  SIGMAS 

COMMON/1/  MSK  (  1000  I  *  KS22  (  80,240,3  )»IV('>  5,4,3)  *X  1  (4,31  ,X2 ( 4 , 3 )  » 
1  X3(4»3)  *  1C 2 1  <  S 2 >  ,KMX(  12*12)  *CDA  (  72 *3 *53  )  . 

COMMON/1/  MSK(IOOO).  KS22 (80*240.3) .1  V(55 ,*.3 ) .XI (4,31 ,X2 ( A .3 » . 

1  X3(4»3 1  •  *21(82)  «KMX  ( 12*12) »CDA (  72*3*53)  » 

2  IN*  MM,  IP,  I  PL »  A21(82)  *  X(32,2)  ,  ICO<4.3,55> 

3  •  CENTO, 3. 55  )  .  OZ 
EQUIVALENCE  ( KS22 ,STRS ) 

DIMENSION  AS22( 60,240) 

EQUIVALENCE  (AS22( 12801 ) .DEL) 

EQUIVALENCE  (S( 19201  I.AS22) 

COMMON  /INPUT/  XNU(2),  E(2)  ,XLAM(2)  ,  6(2)  ,  CMX(6»6,2)  *  DO 
DIMENSION  S(80,80,5> 

EQUIVALENCE  (KS22(  19201 1  ,S) 

DIMENSION  STRS ( 6, 1500 )  ,  DXI12)  ,  DELOOOO) 


PRINT  OUT  DISPLACEMENTS.  6  PAGES 


PRINT  1000 

1000  FORMAT ( 1H1 ,50X . 1 3HDI SPLACEMENTS,// ) 
JCNT  -  0 

DO  15  J-1,102 

JCNT  «  JCNT  ♦  1 

IF( JCNT.LE.18)  GO  TO  14 

PRINT  1000 

JCNT  »  0 

14  JFIR  •  7»(J-1)*1 
JLAST  ■  JFIR  ♦  6 

PRINT  1001.  (K.K-JFIR,  JLAST) 

1001  FORMAT ( 1H  ,7 ( BX , 4HDEL ( • 1 3. 1H)  ) ) 

PRINT  1002.  (DEL (K ) *K-JF IR .JLAST ) 

1002  FORMAT ( 1H  .7 ( 2X . E14. 7 » / ) 

15  CONTINUE 

PRINT  1001.  (K.K-715,720) 

PRINT  1002.  <DEL(K),K-715,720> 

C 

1123  -  0 
DO  406  11*1,9 
DO  406  12-1, IN 
DO  406  13-1,3 
1123  -  1123*1 
C 

C  II  COUNTS  FLOORS 

C  12  COUNTS  BASE  TRIANGLES 

C  13  COUNTS  TETRAHEDRONS  ABOVE  BASE 

C  1123  COUNTS  ALL  OF  THEM 

C 

KZ  ■  0 

DO  405  KJ  -  1.3 
DO  405  KK  -  1,4 
KZ  -  KZ*1 
ISA  -  1 

I  PL  ■  3#(  I V  (  1 2  «KK  « I  3  >  +  (  1 1  —  1 )  »MM-1  )+KJ 


IF( 

IPL.GT.864) 

GO 

TO  305 

GO 

TO  306 

305 

IF  ( 

KJ.EQ.3)  GO 

TO 

307 

I  PL 

■IPL-96 

GO 

TO  306 

307 

0X( 

KZ)>0. 
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GO  TO  405 
306  I  PM ■!  PL 

in  IPL.LE.il  )  GO  TO  300 
I  PL  ■  IPL-10*(!Pl/  96-»l) 

300  CONTINUE 

I  FI  I PM-MSK ( I  PL ) )  409.418.411 
416  IFCKJ.LT. 3)  ISA  ■  1-2*1C0IKK. I J,  12  > 

410  OXIKZ)  »  DEL ( IPL I  *  ISA 
GO  TO  405 

4C9  IFIIPL.EO.il  GO  TO  416 
I  FI  IPL.E0.2)  GO  TO  425 
IMP  -  IPL-1 

412  I F I  I PM-MSK (IMR))  414.419.416 

419  IFIKJ.LT. 3)  ISA  ■  1-2*1 CO  IKK. 13. 12  ) 

415  OXIKZ)  ■  OEL(IMR)  «  ISA 
GO  TO  405 

414  IMR  ■  I MR-1 
GO  TO  412 

416  OXIKZ)  -  DO 

MLD  «  MOO  13*1 IVI I2.KK.I3I-1 ) *K  J  »3*MM ) 

I F I MLD.GT .19.0R.MLD.E0.2)  OXIKZ)  «  0. 

GO  TO  405 

411  IMR  -  IPL  -  1 

417  IFI IPM-MSKI IMR) )  416.415.421 
421  IMR  «  IMR+1 

GO  TO  417 
425  OXIKZ)  *  0. 

405  CONTINUE 

CALL  MXMULTICDAI 1,13,12) .OX.STRSI 1.  1123)  .6.12.1 ) 

406  CONTINUE 

STRS  NOW  CONTAINS  THE  SIX  STRESS  COMPONENTS 
PRINT  OUT  CENTROIDS.  710  PAGES 

JCNT  *  1 
PRINT  1027 
112  »  0 

DO  50  11  *  1.9 
OZD  «  DZ*  <11-11 
DO  50  12  «  1,  IN 
DO  500  13*  1.3 

CENT  I  3  «  I  3  «  12  )  *  CENTO, 13. 12  )  +  DZD 
500  CONTINUE 

112  *  112  *  1 
1123  »  3*1112-1)  ♦  1 
1123P2  *  I123+2 

I 123M1  ■  1123-1 
PRINT  1020 

PRINT  1021.1  I  ll.STRSI l.J) ) » J*I  123.1  123P2 ) 

PRINT  1022,1 ( 1 2. STRS  I  2.  J) )  . J-1 123. I123P2) 

PRINT  1023.  I ( J.STRSI 3.J*I  123M1 ) ) ,J*1  ,3 ) 

PRINT  1024. K  CENT(1,J,I2).STRS(4.I123M1+J)).J=1,3) 
PRINT  1025.(1  CENTI2,J,I2).STRS(5.I123M1*J)),J»1,3) 
PRINT  1025.(1  CENTO,  J,  12).  STRSI6 . I 123M1+J ) ) , J«  1,3) 
PRINT  1026 
JCNT  ■  JCNT* 1 
1F| JCNT. LT. 7)  GO  TO  50 
PRINT  1027 
JCNT  ■  0 
50  CONTINUE 

1020  FORMAT  I 1H  .3 (20X, 13HS TRESS  VECTOR. 6X )  > 
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1021  FORMAT ( 1H  »3 ( 5HLAYEP »9X*I3*3X»E13.6*6X) ) 

1022  FORMAT ( 1H  .3 ( 8HTRIANGLE .6X .  I  3 .3X .E13.6 t6X ) ) 

1023  FORMAT ( 1H  .3 ( 1 1HTETRAHEDR0N *3X » I 3*3X . E13.6 *6X ) ) 

1024  FORMAT ( 1H  .3 ( flHCENTRO I D • 2X t F7.4 • 3X .E 13.6 *6X ) ) 

1025  FORMAT ( 1H  .3 ( 10X.F7.4 . 3X .E13.6 t6X ) ) 

1026  FORMAT (1H  ) 

1027  FORMAT ( 1H1 ) 

RETURN 

ENO 


rtnri  n  n  n  n  n  n  n  n  non  non 


SUBROUTINE  INVERT(B»K»K2.XMIN»FLAG) 


C  THIS  SUBROUTINE  SETS  UP  A  UNIT  MATRIX  ADJACENT  TO  9UiKI 
C  ELEMENTARY  ROW  OPERATIONS  ARE  THEN  PERFORMED  ON  THE  NEW  K  X  2K  MATRIX 

C  TO  REDUCE  B(K.K)  TO  A  UNIT  MATRIX.  THIS  WILL  PLACE  THE  INVERSE  OF 

C  ELEMENTARY  ROW  OPERATIONS  ARE  THEN  PERFORMED  ON  THE  NEW  K  X  2K  MATRIX 

C  TO  REDUCE  b(K.K)  TO  A  UNIT  MATRIX.  THIS  WILL  PLACE  THE  INVERSE  OF 

C  THE  MATRIX  b(K.M  IN  THE  RIGHT  HALF  OF  9(K.2K) 

C  ON  EXIT.  THE  INVERSE  OF  B  RtPLACES  B 

C  B  IS  AN  ARRAY  OF  2*K**2  LOCATIONS  CONTAINING  THE  MATRIX 
C  K  IS  THE  DIMENSION  OF  THE  SQUARE  MATRIX  B 
C  K 2  IS  2*K 

C  XMIN  IS  THE  SMALLEST  ALLOWABLE  MAGNITUDE  OF  THE  PIVOT 
C  FLAG  WILL  BE  RETURNED  AS  0.  IF  THE  INVERSION  WENT  OFF  OK 
C  FLAG  WILL  BE  RETURNED  AS  1U.  IF  A  PIVOT  ELEMENT  WAS  TOO  SMALL 
C  FLAG  SHOULD  BE  TESTED  AFTER  EACH  CALL  TO  THIS  ROUTINE 
C 

DIMENSION  B ( K.K2  ) 

C 

FLAG  *  0. 

SET  UP  UNIT  MATRIX 

DO  I  1=1. X 
DO  1  J=1.K 

BCI. K+J)  =  0. 

I F ( I.EQ. J)  B ( I ,K+J )  =  1. 

1  CONTINUE 

FIND  LEADING  ELEMENT  WITH  GREATEST  MAGNITUDE 

DO  6  J=1.K 

M  =  J 
N  *  J+l 

IF(N.GT.K)  GO  TO  21 
DO  2  L=N.K 

IF  ( ABS IB(M.J)  )<LT.ABS(B(L.J) I )  M=L 

2  CONTINUE 
21  CONTINUE 

IF  ( ABS (B(M.J) 1.LT.XMIN1  GO  TO  10 

INTERCHANGE  JTH  AND  MTH  ROWS 

DO  3  L=J»K2 
D  a  BIJ.Ll 

BCJ. L)  =  B(M.L) 

B(M.L)  =  D 

3  CONTINUE 

ZERO  OUT  PIVOTAL  JTH  COLUMN.  SKIPPING  PIVOTAL  JTH  ELEMENT 

DIVIDE  JTH  ROW  BY  PIVOT 

DO  4  M=N.K2 
B(J,M)  =  B(J.M)  /  B(J.J) 

4  CONTINUE 
DO  6  M-l.K 

M  DETERMINES  ROW  BEING  MODIFIED,  ONE  WHOLE  ROW  AT  A  TIME 

IF  (  M.EO.J  >  GO  TO  6 
DO  5  L*N.K2 
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c 

C  L  DETERMINES  ELEMENT  IN  THE  MTH  ROW 
C 

B(M.L)  «  8(M*L )  -  B(M*J)  *  B(J.L) 

5  CONTINUE 

6  CONTINUE 
C 

C  INVERSE  OF  B  IS  NOW  IN  RIGHT  HALF  OF  BU.K2) 
C  NOW  MOVE  B  INVERSE  TO  WHERE  B  WAS 
DO  7  I  *  1  *K 
DO  7  J*1.K 

B (  I « J  )  *  B ( I  * J+K  ) 

7  CONTINUE 
RETURN 

10  FLAG  *  10. 

RETURN 

END 
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SUBROUTINE  MXMUL 1 ( A .B .C ,M ,N ,K ,M1  »N1  .LI) 

C 

C  THIS  SUBROUTINE  MULTIPLIES  MATRIX  A  BY  MATRIX  B  AND  STORES  THE 
C  PRODUCT  IN  C.  (C  CANNOT  BE  THE  SAME  AS  A  OR  B. ) 

C  A  IS  (M  X  N) 

C  B  IS  IN  X  Kl 
C  C  IS  (M  X  K) 

C 

C  Ml  IS  NUMBER  OF  ROWS  !N  ARRAY  A  IN  CALLING  PROGRAM,  Ml.GE.M 
C  N1  IS  NUMBER  OF  ROWS  IN  ARRAY  B  IN  CALLING  PROGRAM,  Nl.GE.N 
C  LI  IS  NUMBER  OF  ROWS  IN  ARRAY  C  IN  CAlLING  PROGRAM,  Ll.GE.M 
C 

DIMENSION  A (Ml »N )  ,  BiNl.K)  ,  C(L1  ,K ) 

C 

DO  1  I«1»M 

DO  1  L*1.K 

C  {  I  » L )  =  0. 

DO  1  J* 1  ,N 

C ( I »L  )  *  C(I.L)  +  A ( I , J )  *  B ( J  »L  ) 

1  CONTINUE 
RETURN 
END 
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SUBROUTINE  MXCON(A.B.X«M.N> 


C 

C  THIS  SUBROUTINE  MULTIPLIES  MATRIX  A  IMXN)  BY  CONSTANT  X.  RESULT  IN  B 
C  A  MAY  BE  SAME  AS  B. 

C  THIS  SUBROUTINE  MULTIPLIES  MATRIX  A  (MXN )  BY  CONSTANT  X»  RESULT  IN  B 
C  A  MAY  BE  SAME  AS  B. 

DIMENSION  A(M.N)  •  B(M.N) 

00  1  I  «1 *M 
DO  1  J«1.N 
B( I t J) *  X*A ( I .J) 

1  CONTINUE 
RETURN 
END 


nonnnnnn 
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SUBROUTINF  MXMUL  T  (  A  .  B  .C  tM .N .K ) 

THIS  SUBROUTINE  MULTIPLIES  MATRIX  A  BY  MATRIX  B  AND  STORES  THE 
PRODUCT  IN  C.  (C  CANNOT  BE  THE  SAME  AS  A  OR  B. ) 

A  IS  (M  X  N) 

B  IS  (N  X  K ) 

C  IS  (M  X  K) 

DIMENSION  A(M.N)  •  8 ( N  *K I  •  C(M.X) 

DO  1  I  =  1  *  M 
DO  1  L* 1 *K 
C(I.L)  *  0. 

DO  1  J*1.N 

C(I.L)  »  C(I.L)  +  A ( I  * J)  *  B( J«L ) 

1  CONTINUE 
RETURN 
END 


nnnn^nnn^ 
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SU3R0UT INE  MXSUB ( A, B.C *M ,N ) 

THIS  SUBROUTINE  SUBTRACTS  MATRIX  B  FROM  MATRIX  A. STORES  RESULT  IN  C 

THIS  SUBROUTINE  SUBTRACTS  MATRIX  B  FROM  MATRIX  A, STORES  RESULT  IN  C 

THIS  SUBROUTINE  SUBTRACTS  MATRIX  B  FROM  MATRIX  A. STORES  RESULT  IN  C 

A*  B.  AND  C  ARE  (M  X  N)  (C  CAN  BE  THE  SAME  AS  A  OR  B) 

DIMENSION  A(MiN)  .  B(M.N)  .C(M.N) 

DO  1  I-l.M 
DO  1  J *  1  *N 

C(I.J)  «  A ( I • J )  -  BII.JI 
1  CONTINUE 
RETURN 
END 

1056  CARDS 


COMPUTER  PROGRAM  FOR  POINT -MATCHING  METHOD 


SRp4  SnEFF  ~  - - - - 

ST/  J5T$P/_59T$C/200ES 

SFTN.L.P.'  - -  - - 

_ PROGRAM  PMATCH 

C  THIS  PROGRAM  READS  THE  INPUTS  AND  SERVES  AS  A  DRIVER  FOR  THE  SUBROUTINES 
C  SOLVING  THE  TWO-DIMENSIONAL  STRESS  PROBLEM  OF  A  FIBER-REINFORCED' SAMPLE 
C  USING  THE  POINT  MATCHING  METHOD.  THE  EQUATION  SOLVED  IS  C.X  =  Y  ,  WHERE 
C  G  IS  INEQ'  X  NUNK )  »  X  IS  (NuNK  X  1).  AND  Y  I S  <  NUNK  'X  IT 

C  SHOULD  M  AND/OR  MPR  BE  CHANGED.  DIMENSIONS  OF  COEF  SHOULD  BE  8E-ESTABL I SHFD . 
C  AS  WELL  AS  DIMENSIONS  OF  6"  AND  ON .  AND  NEuVnX.AND  NUNK 

_ COMMON  / 1  /  _G£1_48  ,S7  )  .  GNI56.57)  ,  E(2>  .  NU  ( 2  )  . 

l  a'i'n.  B'lN.  CIN.  VF 

_ DIMENSION  ERR ( 148 )  .  Y(S6)  .  GNSAVIS6.S7) 

Type  rea'l'nu -  - - - - - 

_ COMMON  /SIZES/  NEC,  m,  ,v,pp,  nx  ,  NUNK  .  NUNKP1 

DIMENSION  COMENT 1101  . 

_ DATA  ( PI_=3  ._141592  7 )  t  (NEO=148)  .  (M  =  7)  ,  (MPR=7>  .  (NX  =  37) 

DAT  ATE  LM  IN=1 • E-8 )  .  (NUNK=56),  (NUNKP1=57) 

_ JDUM  =  IN1T< 1.1J 

C  IN  I T  IS  CALLED  HERE  TO  INITIALI7E  BINOMIAL  COEFFICIENTS  TO  ZERO 
_ KASENO  =1 

READ  1 000 .(COMENTIi). 1=1.10)  - 

1000  FORMAT ( 17A8 ) 

r'RFA^lTrOTiETTlV  ET2 J  .  NUflt,  NUI2).  VF,  BIN  - - - - 

1001  FORMAT  ( 6F1 0,4 )_ 

IFIBIN.EO.O. )"  GO  TO  2  '  "  " 

_ PRINT  2002  .  K_A S ENO 

2002  FORMAT ( 1 3H1  CASE  NUMBER .12) 

£Rffl,O200(!Affc$MENl(  I  >  ,1  =  1,10)  . . . — . . 

2000  FORMAT tlxluAfl  )_  _ 

CIN=SORT( 3. ) /2.*BIN  "  -  - “ - 

_ AIN-  SORT ( 4. »  CIN*BIN*VF/PI ) 

PRINT  200  1  ,  "e"(  1),  ~E  (  2  )  .  Null').  NU(2).  VF ,  AIN.  BIN,  CIN 

2001  FORMAT (7X  2HEI  , 1 2X ,3HF I  I » 1 IX ,4HNU I  .10X .4HNUI I .11X.2HVF.13X, 1HA , 

1  13X.lHB.13X.iHC  /1X.8I2XE12.S)  )  - - 

_ _ _ 

C  NEO  -  NUMBER  OF  EQUATIONS  IN  OVER-DEFINED  SYSTEM 
C  M  -  SUMMATION  LIMIT.  REGION  II 

C  MPR  -  SUMMATION” LIMIT.  REGION  I  - - - - 

JC _ NX-  TOTAL  NUMBER  OF  UNKNOWNS  ,  REGION  I 

C  NUNK  —  NOmBER" OF  UNKNOWNS  RFGION  I  PLUS  REGION  II.  . . . 

C  NUNKP1  -  NUMBER  OF  UNKNOWNS  PLUS  1 

c  Elm i n  -  mImmum  allowable  magnitude  for  a  pivotal' 'ElEmEnt - 

_c _ _ 

C  '  - - 

C  SET  UP  THE  MATRIX  G ( NEUXNUNK ) ♦  THE  COEFFICIENT  MATRIX 

C  IN  THE  OVER-DEFINED  SYSTEM  GX  =  Y  .  AND  ALSO  SET  UP  Y,  THE  INDEPENDENT 
C  VFCTOR  _ 

CALL  GETG  —  -  - 

_C _ _ 

C  NOW-  HAVE  KNOWN  MATRIX  G  AND  INDEPENDENT  VECTOR  Y.  GET  PRODUCT"  GN  = - 

E.  IG-TRANSPOSF)  *  G.  GN  IS  (NUNK  X  NUNK)  . 

c  ~  .  . . 

_ CALL  MXTMUL(G.G.GN.NEQ.NUNK.NUNKP1 ) 

C _ NOW  HAVE  LEAST-SQUARES  COEFFICIENTS  IN  GN  (NUNK  X  NUNKP1 ) 

C  "  -  -  -  - - 

C  SAVE  GN  FOR  BACK  SUBSTITUTION 
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c 

_ J30  9100  I  =  l  ,NUNK 

DO  9100  J=1 .NUNKP1 
9100  GNSAV(I.J)  =  ONII.J) 

C 

_ CALL_F0NSLV(GN,NUNK.NUNKP1 .FLAG.PLMIN) 

"  Tf  (FLAG.NF.O. )  GO  TO  100 


x _ _ 

C  MOW  NlJNKPl  TH  COLUMN  CF  GN  CONTAINS  UNKNOWN  VECTOR 

~c  caUculaTe  ERRORS  I N  primitive  system  and  "print  them 


PRINT  6070 

6070  FORMAT  (__  1_H1 , 2GX . 26HE RRORS  IN  PRIMITIVE  SYSTEM///) 

"on  Iso  T-T.'iTb 

CALL  VECMUL  1  G»GN  ( 1  .NIINKDl )  .ERP(  I  )  .NEO.NUNK._U 
ISO  FPR(i)  =  GCI.NUNKP1)  -  FRRII) 

DO  200  I  =  1.21 

IF  IP  =7*(T-l7*i 

I  L  A  S  T  =  6  +  IFIR 

. —IFH.E0.i5r  PRINT  6070  ' 

PRINT  607?.! J.J=IFIR,ILAST) 

200  PRINT  6071, ( ERR ( J) . J= IF IR , I  LAST  ) 

PRINT  6077  .  IAS 
'dpTNT'6'07'1  ,  ERR  (T4'8  ) 

6071  FORMAT! 1H  .7I3XE13.6)/  ) 

"&0T2  "FORWrriH-V7TIOX  ."ZFOCTVIT.lfn  5 )  '  )  ~  . 


C 

r 

c 


'SYSTEM' 


ppi'nt'  60R‘o- 

6080  FORMAT! 1H1 .2UX.2AHERRORS  IN  SQUARED  SYSTEM///) 

- .  . . 


CALL  VECMULIGNSAV.  GN! 1 .NUNKP1 ) .ERR ( I ) .NUNK.NUNK.I ) 

- lSCERRTIT  =  'GN'SAV!  I  .NuNKPTT -ERRTn - : - 

DO  190  1=1,8 

if7r  =  7*(f-rr+r  . 

I  LAST  =  1FIR+6 

. 'PR'l  N T" "  6 07  2Vf  J  rj= TFTR  ,T L  AST  j 

1O0  PRINT  6071. (ERR! J).J=  IFIR.ILAST) 

c  substitute  unknowns  into  equations  f o r' s trEssES  and  Displacements. 


c 

.  call ‘ba"cksb'“  . . 

c 

~  C  XlCSTRESsES'  And'  D'l SPlaCEMEntS  are  known.  GET'Mo6uluS''EE~/~~ETT 


GO  TO  2 

100  PPINT  200S,  FLMIN 

2 0 O's'TORM AfTSfeH  1A"L A R G E  S  T  PIVOTAL  ELEMENT  WAS  SMALLER  1  N'^GRTf T U  6  P ' YH AN 
1  .Eli. 4) 

2  cTonTTnUF 

STOP 

END  . .  . . 


/ 
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SURROUT  INF  GETC, 

COMMON  /l/  0(148.571  .  ON ( 56 . 57  ) 

1  aTnV  bin'.'cTnT  jf 

TYPE  REAL  NU  _ _ _ _ 

COMMON  /SIZES/  NEC.  M .  MPR .  NX.  Nl 


THIS  SUBROUTINE  GENERATES  THE  COEFFICIENT  MATRIX  (0  MATRIX) 
OF  THE  KNOWN  OVERDEFINED  SYSTEM 


TYPE  REAL  K  A  Y _ 

DATA  <PI=3. 14159271 
KAY  «-E( 1  )/( 1«+NU(  1  )  I 
ZER*0  OUT  'AUGMENTED  MATR1 
_  D0_1_I±1  .NEO 

do  i  j*  iTnunkp  l” 

r.  1 1 1  j  ) «  n  ■ _ 

1  CONTINUE 
EQUATIONS  -)r‘  TO  35 

x x'*  cTn'-  aTn  . 

YY*BlN/2. 

“C  SU  L*  S'l  Gift  Y.  T»  X  X  VY  Y  .  1 
CALL  SIGXI 1.2.XX.YY.- 
CALL  TAUXY(2.1.XX.YYi 
CALL  TAUXYI3.2.XX.YY, 

call"  (V'.T.VxTyV;  i  ,'T 

CALL  U(4,2ji_XX.YY.-1.  I 

'TA"ll"v  T5VI  .XX.  V  y  .1 .  V 

CALL  V(6,2,XX.YY,1.  ) 


EQUATION  ^6 

DO  1 0  *J=1 - 14 

__cpj  *cn  sj_p  I  *  j  /  3  u  - 1 
SPj“sfMPT»'J7To.  i 

SPJ2  =  5IN(°I*J/15.  )  __ 

XX  =  CIN-AI\'«CPJ 
Y Y  =  B1 N A2  : N*SPJ 

CAL  L  s'i  gx  (  J+'&7 17  X  x'.VyVc  P  J  *  *  2  > 
__  CALL  SIGXIJ+6.2.XX.YY . — CP J**  2 ) 
CAL L~S  I G y'(  J+6 l7 XX  ,  Y Y.SPJ**?  ) 
CALL  SIGYI J*6.2.XX.YY,-SPJ««2) 
call  TAUXYI J+6.1  .XX.YV.SPJ2) 

_ CALL  TAUXYI J+6.2.XX.VV.-SPJ2  ) 

Tquat'ion  37  ~ 

CPJ2=  C0S(PI*J/15. ) 

SPJ2  =  SPj2/27 . 

call  SIGYIJ+2U.1  .XX.YY.SPJ2) 
CALL  SI  GY  I J+20.2.XX.YY.-SPJ2 ) 

_ CALL  SIGXI J+20.1 .XX .YY.-SPJ2 ) 

'""call  "sTg xT j + 2 o'. 2 »  XX  *  Y yV S  P  J  2  ) 

_  CALL, TAUXYCJ+20. 1 .XX ,YY .CPJ? ) 
"CALL  "TAUXY  (  J+2  h  ,  2  .XX,  Y Y  ,  - C P  JZ  ) 


EQUATION  J8 


CALL  U< J+34,1 .XX.YY.CPJ) 
call  IJ<  J+34.2.XX.  YY.-CPJ) 
'c'al'l'V  "f j+34,  i  .  XX  ,  YY  »  SPJ ) 
CALL  V( J+34.2.XX.YY.-SPJ) 


EQUATION  39 


CALL  IK  J+48.1.XX.YY.SPJ) 
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call  U( J+48*2.XX.YY.-SPJ) 

_ CALL  V(J*48»1. XX.YY.-CPJ) 

CALL  V<  J*4fl,2,XX.YY.CPJ> 
1A  CONTINUE 


E0U4TI0NS  67  TO  68 


XX-CIN _ 

’Y"Yi"BIN"7yr-Ai'N" 


» 

CALL  SIGY<67.?,XX.YY,-1.  ) 
CALL  TAUXYI64, I ,XX, YY.l  . ) 

f 

CALL  TAUXY(65.?.XX.YY.l.  ) 
CALL  IH66.1 .XX. YY.l  .  ) 

f. 

CALL  UI67.2.XX.YY.1.) 

CALL  VI68.1.XX.YY.1.) 

i 

CALL  V(6B.2.XX,YY.-1.  > 

G ( 66  .NUNKP 1 )  *  KAY 

l 

c 

0(67. NUNKP 1  )  «  KAY 

... 

ToUATfoNS  46  7CnI>  a  7 

Bo  20  J-1.9 - 

XX=CIN-AIN+J*AIN/10. 

TV  .077575. 

CALL  TAUXY(J+66.1.XX,YY,1.) 

c 

CALL  v ( J+77.1 .XX.YY.l.) 

c 

23  CONTINUE 

EQUATIONS  48  TO  53 

c 

XX-CIN 

YY-BlN/2. 

CALL  TAUXY(B7*l*XX*YY»l»i 

Ct Al L  U ( §8  *  1  »XX  * YY  *  1 •  t 

CALL  VIBO.l.XX.YY.l.) 

c 

G ( 8 8 i NUNKP 1 )  «  KAY 

EQUATIONS  51  AND  52 

C 

DO  25  J-1.9 

VY-BIN/2.-JMAIN/10.) 
call  TAUXY(J+89,1,XX.YY,1.) 

r 

CALL  U< J*98,l.XX,YY.l.) 

25 

CONTINUE 

£ 

C 

EQUATIONS  53  AND  54 

DO  50  J-1.4 

YY»BIN/2.-AIN-J*(BIN-AIN)/5. 

" call  TAUxviJ+r^.yrsre.YY.i.i 

CALL  U( J+lll .2.XX.YY.1.) 

33 

Gl J+l 1 1 .NUNkPI )  -  KAY - 

CONTINUE 

c 

T 


EQUATIONS  55  TO  57 


YY--B1N/2. 

CALL  TAUXYl  116.2  tfcx.'YY.l  .  f 
_CALL  U(  U7. 2, XX. YY.l.  > 

gTTT7  .nunkPTT”"-- kA"y 

CALL  V(  118. 2. XX. YY.l. ) 
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c 

C  EQUATIONS  58  AND  59 

C  . 

DO  AO  1-1.12 _ 

XX  ■<  I-1I  •ICIN/12.  ) 

YY«-B I N/2 « 

c a'l’l' t‘aux yTi+  IY8T2 Vx x . yy . f . ) 

_ CALL  V ( I »1 30 1 2.XX.YY.1 .) 

40  'continue'" 

EQUATIONS  60  AND  61 

5o'"4TT*t;v  t 

XX  «(  I-n  *£C IN-AIN)  A3. 

*YY«BlN/2.  '  . 

CALL  TAUXYI I+142.2.XX.YY.1.  ) 
CALL  VI 1+145 .2. XX. YY. 1 . ) 

45  CONTINUE 


C 

C 


t 


t 

I 

! 
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SUBROOT  INF  BACKSB 

COMMON  / 1 /  G  (  1 4  8  » S  7  I  ,  GNI56.57)  ,  E(2)  .  NU(2>  ♦ 

1  AIN.  BIN.  CIN,  VF 

_ TYPE  REAL  NO  _  _ _ 

COMMON  /SIZES/  ~NEQ,  M.  MPR.  NX.  NUNK  ♦  NUNKP1 
DIMENSION  DF LP ( 2  )  .  AyE(2)  .  SIGNI20.2I  .  TAONI20.2) 

DIMENSION  JI (11.11) 

DIMENSION  XFI2o.ll)  .  YFI20.il)  .  STRESS ( 20. 1 1 . 5  )  .PHI ( 20 ) 
DIME NS  ION"  PH  I  DEG ( 2  u  ) 

_ TYPE  REAL  Jl  ' _ 

C 

TYPE  REAL  KAY 

. . "DATA  “T(  AYE  i  I  )  .  I  =1 .2  )  =2H  I.2HII) 

KAY  =  -F  (  1)  /(  1  •  +  NU  (  1  I  ) 

"c"""deve"lop""valu"es  or  sigmax,  SIGMAY  on  10X10  GRID,  STARTING  at"  x=o. 

C  Y  «  B/2.  .  WALKING  along  ROrf  TO  RIGHT.  _ 

c 

DtLX  =  Cl N/l w 

. . DEL  Y  =”-B  IN  /  10.  " 

XI  =  -DELX 

fppN-f'  =  '0‘ 

DO  10  J  =  l,ll  _ _ _ _ _ 

I P  R  N  T  =  IPRNT  +  1 
GO  TO  (21.22 )  .  IPRNT 

2T~'p rTnt"' l'67f i"~  . . 

GO  TO  23 

'""T2"PRWT0T0 -  - - - . . 

IPRNT  =  0 

23C6NY1NUE . 

XI  =  Xl+DELX 

Y-o“;_bTn/2  I-WlY 

DO  10  1=1,11 
. 75  ="yo"+"df"l  Y 

XF ( I , J )  =  XI 
YF ( I  ,  J )  =  YO 

I F  (  KIN-XFI  I  ,_J)  )**2+(_BIN/2.-YF  I  I  ,  J)  )  *  *2  .  GT  .  A I  N  *«  2_) 

Too  "to” 

L  =_  1 

-j-j-jj-j,  -  a‘yeTI")  " 

_ GO  TO  6 _ _ _ _ _ 

5  L  =  2 

_ J  I  ( I « J )  =  AYE ( 2  ) _ 

-  contTnTje  . 

C 

C  K‘f  TELLS  WHK H"HA L r OF'-UNKNOWN  VECTOR  TO  USE"  "IN  MULf  TPITl CATION" 

C _ 

C  ZERO  OUT  MATRIX  G  AS  REQUIRED 
C  _ 

DO  7"l2“l  =  1.5 
DO  7  J21  *  l.NUNK 

7  g  j“f2 1",  J  217"  =  "  0"T  "  . . 

C _ _ _  _ _ 

CALL  SIGX(1,L»XF(I,J),YF(I,J),I.) 

CALL  S I  GY ( 2  »  L , XF ( I , J ),YF ( I . J )  ,  1 .  ) 

C"AL’L"‘fA"uXY"r3'ir,x'F'('fVj(  ^YF('i  rj)  ,  l".  ) 

CALL  U(4,L.XF(I,J),YF(I,J)«-1./KAY) 

c"al"l-"v1'5  7l  Vx’FTr.'jTiYFTr,  J")",  - 1 .7*  ay'  i 

c _  _  _  _ _ 

DO  9  IM=  1 , 5 

9  CALL  VECMULI G.GN ( 1 .NUNKPl ) ,  STRESS ( I ,J , IM  )  ,  NEO.NUNK.IM) 

c"  pri  nT'ouT "  g'rTd""v a'rTab  l’e  s 

c 
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"p ri’nT'T 6  o  i  ,  x  f  i  i  ,  o  i  i  vc <  i  *  j  i  .  i strfssi  i  .o.kT.km  rsTrjTH  ijj 

in  CONTINUE 

c 

C  STRESS ( I  * J  » 1 )  CONTAINS  SIGMA  X  ,  1  =  1.11  .  J  =  1  *  1 1 
C  STRESS! I »0.2T  CONTAINS  SIGMA  Y  ,  I  =  l.]l  .  J  = 1 . 1 1 

CL  _ET„CETFR_A _ 

C 

1000  FORMAT  C  1H1,13X,1HX,13X.  IHY.fiX.AHS  IGM  AX  ,  BX  «f>HS  I GMAY.9X.5HTAUXY.13_X. 
F"  "  lHU.nX.  1HV.  10X.6HREGION/ ) 

1001  FORMAT  (  1H  2I3X.FU.4)  . SOX. Ell. A)  .12X.A2  /)  .  _ _ 

1070  FORMAT!////) 

C 

C  now"  MOVE  ALONG  EDGE  OF  FIBER. FROM  THE  TOP.  COUNTER-CLOCKWISE. 

C _ 

DELPHI  =6. 

_ PH  1 1  =  0. _  _  _ _ _ 

DO  0  7  1  =  1,14 
r>HI  1=dhI  l+OELPHI 
PHIOEG  (  I  f  "=  PHI  1 
PHI ( I )  =  PHI  1  /57. 29578 
.  x'FTi  .1  )  =  C  IN-AIN*C0S  (  PH  I  (  I)  ) 

_ Y F ( 1,1 ) =  B IN/ 2 . -A  I N*S I N ( PH  I  ( I  )  )  _ 

.00  29  J=1  .2 
C.  . . 

C  YfRO“oI)T  MATRIX  G  AS  RFOU 1  RED 

c 

DO  28  121*1,5 

_ DO  28  J1  =  1  .NUNK  _  _ _ _ 

28  G!  121,01)=' 

CALL  SIGXU.  J,  XK  I  ,1  )  ,YF(  I  ,1  )  ,  1 .  1 

C  ALiT  S  IGY  (2.  J.  XF(i,l),YFf  l.lt.l.  ) 
call  TAUX Y 1 1 i Ji  XF( I ,1 )  ,YF(  I ,1  ) ,1 .  ) 

CALL  U ( 4. J.  XF(  I  ,1 )  ,YF<  I  ,1  )  ,-1,/KAY) 

_ CALL  V!S.O.  XF(  I  ,1  )  ,YF  !  I  .1  )  ,-1,/KAY’  _______ 

DO  29  K=1 ,5 

29  CALL  VECMUL  ( G.GN ( 1 .NUNKP 1 ) ,  STRESS!  I  »K,J> .NEOjNUNK ,K  ) 


< 


I 


C. 

r 


c 

SIGMAX 

J 

IN 

STRESS! I ,1 ,J ) 

f 

0  =  1.2 

c 

S IGM ay 

J 

IN 

STRESS! I ,2  » J ) 

♦ 

J=1  ,2 

c 

TAIJXY 

J 

I  N 

STRESS! I.3.J) 

* 

J=1  ,2 

c 

U 

J 

IN 

STRESS! I ,4 , J ) 

* 

0  =  1,2 

c 

V 

J 

IN 

STRESS!  154, J) 

* 

0=1,2 

‘c 

c 


DO  SO  J  =  1  ,2 

SIGN! I . J)  =  STRESS!  1 . 1 . J )*COS! PHI ( I ) ) **2  + STRESS (  I ,2,J) 
-  *  SI N( PHI ( I )  )**2  +  STRESS< I  ,3,JI*2. 

_?  *  SIN!  PH  1(1))*  COS!  PHI  (  I  )  ) 

30  TauniT.—;  (STRESS!  I  ,2  ,  J  ) -STRESS  I  I  .1  ,J)  )*SIN(PHI  (  I  )  ) 

_ 1  _  *  COS ( PH  I ( I )  )  +  STRESS!  I  .3. J) 

2  . . .  “*(COS( PH!  (  1 ) ) **2  -  SIN ( Ptll  ( 1 1 ) **2 ) 

33_  CONTINUE 

"c  . . 

c  NOW  HA VF  values  ON  I  NT EPF ACF ,PR I N T  THFM  OUT 

c  ' 

100S  FORMAT! 1 H 1 ) 

PRINT  1037 

1037  FORMAT! 1H1.1SX.19HVALUES  ON  INTERFACE,////) 
b7)"44~J=l  ,2 
DO  44  1=1.14 
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IF(l.FQ.l)  PRINT  1010 

_ PgXN  T  10JL1  j  _  A  YE  (  J  ) ,  PH  1  DEG  (  I  )  ♦  XF  <1  .  \  )  ,  YF  (  I  ,  1  )  ..SJ.GN  ( I  *  J )  *  TAUN  {.I  .  JJ 

1  .STRESS!  I  *4  *  j)  •  STRESS!I,5,J) 

1010  FORMAT ( 8H  REGION  ,  8X.3HPHI  ,  1 3  X  * 1 HX . 1 4 X 1HY ,  1 1 X  ,7HS IGMA  N,  9X, _ 

1  5HTAU  N.11X.1HU.14X.1HV/)  . . 

IFII.no. 14)  PRINT  1012 

Ic 12  F OR MAT <7 77  71 

44,  CONTINUE 

1011  "FORMAT  (Yh”"4X  » A2  » 4X»F10.2*»5X ,6 1 3X  »E12. 5  )  ) 

■  ■  C _ _ _ 

C  NOW  FIND  VALUE  OF  P 

_ £. _  _ _ _  _ 

C  OELP(J)  —  INTEGRATION  STEP  SIZE  IN  REGION  J 

_ .t _  _ _ _ 

P  =  0. 

_ El  =  _ _ _ _ _ 

DELP(l)  =  AIN  /  20. 

CELPI?)  *  (fl IN-AIN)  /  20. 

yY  V 'rTn72*.“dTlpTi1 . 

DO  60  J*I  ,2 

----- . 

C  J«l.  REGION  I  .  J  =  2,  REGION  II 
C 

DO  60  I  =  1,20 

~c . . . .  . 

r  I  STEPS  ALONGJLINE  x=c 

i F ( I.EQ.1.AND.J.E0.2)  GO  TO  60 

- yr  V-VT  _  m--Lp(j) .  . .  - . . . 

C 

- r - - - - 

c 

- OT--5D--j5T-VI  vNuiqr -  - - - 

50  G( 1 , J?  1  )  *  o. 

- raur'57  6x  ( IVJ.'OWVVT.IV)  •  - - 

CALL  VECMUL (G«GN(1»NUNKP1),SX»NEQ»NUNK»1 ) 

c' .  . .  ~ 

C  INTEGRATE  SX  FOR  P.  USING  SIMPSONS  RULE. 

— -f--n"HFTNTK'RTSrT^“3rrT6uAL'T0'"ZERO  ON  'THE  "FTPr5T-PAT5V - 

C  TRAPEZOIDAL  INTEGRATION  IS  USED  ON  THE  SECOND  PASS 

- C  "whEn" 'JTTdrrT'RbM  ' "ONE"  TO  ‘twO,  Th>  I NTE'G RTTf  T(5N— STEP-  "CRXNG'ES  "50  THAT - 

C  AGAIN  SIMPSON, S  RULE  COMMENCES  ON  1=3. 

"  c  .  . 

IFII-2)  56,52,53 

- T?-T'V~pr*T;'13rCFrrrr*T5X+TXTT72.  - - - 

GO  TO  55 

- 53 ""P  P2' '  +  DELPTJ I-4TSX2+'  4.*SX1+  SX")  7~ T. - 

55  CONTINUE 

c  . ' .  .  . 

C  SAVE  PAST  VALUES  OF  P,SX 

'c .  . . 

P2  =  PI 

- pl  p - 

56  SX2  =  SX1 

s"xT“§x' . . . . .  . 

PRINT  9091,1 ,P. PI, P2 

9-09Y~F3SmXT1"Sh  "1V;T3V5hP~=  »"El3.6 »  3hP1  ■ .E13.6 ,5WP7iy£  Y376‘) 

PRINT  909?,  SX,SX1,SX? 

909?  FORMAf (4H  SX  =  , El  3 . 6 ,4HSXl=  ,T\3 .6 , 4HSX 2= ,El3-6  ) 

60  CONTINUE 

--  ~c~  p  _  -p-R-rslNf  value  or  Integral 

C  PI  -  1ST  PAST  VALUE  OF  P  _ 
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C  P2'--''2ND"PAS'T  VALUt'OF  P 
C  SX  -  PRESENT  VALUE  OF  SIGMAX 
C  S_X  1  — ■  l"S T  P/Tsf* VALUE  OF"  SIGMAX 
r  SX?  -  2ND  PA5T  VALUE  OF  SIGMAX 
C 

C_ _ 

C  NOW  HAVE  P.  PRINT  IT 

C 

PRi'NT  102U.P 

in?0  FORMAT ( // / /4H  P=  tE15.8) 
RETURN 
END 


nolo  |  njn  nlo  nj  !  nlon!  |  !  nnr.|n  |.-i  n  run  olo  n  r>  o;n  r>jO  o.o  j  J  .  o;nor,o.r, 


t 


subroutine:  vecmul(a.b,c.i»j«<> 


This  subroutine  multiplies  the  <th  row  of  matrix  a  times 

_ THF_COLII«M  VECTOR  STORED  IN  B  AND  STORES  THE  RESULT  IN  C. 

MATPIX  A  HAS  I  ROWS 
"matrix  "A  HAS  J  COLUMNS 

^  DIMENSION  A(  I.JI  ,  B(J> 

*c  =  0". 

DO  1  L  =  1,J  .  _ _ _ _ 


1  C  *  C+A(<.Li*S(Ll  _  ..... _ 

RETURN 

E  NO 

. SUBROUTINE'  EQNSLVI  B  ,K  »KP1,  FLAG, ELMIN) 

This  subroutine  does  a  gajssian  elimination  procedure  on  the 

AUGMENTED  MATRIX  B  =  AY  ,  WHFRE  AX  =  Y.  X  UNKNOWN. 

THE  INDEPENDENT  VECTOR  Y  STORED  AS  THE  K+1ITH  COLUMN  OF  R  • 

P  IS  (KXK+l).  THE  SOLUTION  ALGORITHM  PROGRESSES  ACROSS 
"Columns  to  the  RIGHT,  CHOOSING  THE  COLUMN  ELEMENT  WITH  LARGEST 

MODULUS  AS  the  PIVOT.  IF  THIS  ELEMENT  IS  LESS  THAN  ELMIN,  FLAG _ 

Ts— STT  TO  1  n .  AND  A  RETURN  TO  THE  CALLING  PROGRAM  IS  EXECUTED. 
FLAG  IS  RETURNED  AS  0.  IF  NORMAL  COMPLETION  OCCURS.  FLAG  SHOULD 
ALWAYS  BF  tested  ON  RETURN  TO  CALLING  ROUTINE. 

KOI  IS  K+l.  REQUIRED  FOR  VARIABLE  DIMENSIONING 
"ON  RETURN.  SOLUTION  VECTOR  IS  IN  LAST  COLUMN  OF  B  ,  REPLACING"?". 


)  IMF NS  I OM  RTK , K o 1 f 


FLAG  =  0. 


FIND  LEADING  ELEMENT  WITH  GREATEST  MAGNITUDE 


DO  6  J  =  1«K 
M  «  J 

N  =_J+1 _ 

DCf  2  L*N  »K 

IF  (ABS(B(M«J) ).LT.ABSIB(L,J) ) )  M«L 
2  CONTINUE 

IF  (ABSIBIIM,  J)J  .LT.  ELMIN)  GO  TO  10 
INTERCHANGE  JTH  AND  MTH  ROWS 


DO  3  L=J,KP1 
C  =  B(J,L) 
B(J.L)  =  B(M,L) 
Y<M.L7"'=~D 
3  CONTINUE 


ZERO  OUT  PIVOTAL  JTH  COLUMN,  SKIPPING  PIVOTAL  JTH  ELEMENT 
DIVIDE  JTH  ROW  BY  PIVOT 


DO  A  M=N,KP1 

B7  JTmT“="B'(  J,M)  /  B(J.J) 
4  CONTINUE 
Uo~b  M5T,"< 


DETERMINES  "RO_W”BE  I  NG  MOD  I F  I  ED  ,  ONc  WHOLE  ROW' At  A  TIME 


MTH  ROW 
*  Bl J.L  ) 


IF  (  M.EO.J  )  C,0  TO  ft 

DO  ft  L=N.KP1 


C  L  OF  TFRN I  NFS  FLEWFNT  IN  THF 
B(M.L)  =  B (M . L )  -  B(M.J) 

_ ft  CONTINUE 

ft  CONTINUE 


RETURN 

n  FLAG  =..1Q... 
RETURN 
END 
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FUNCTION  COFF { I ,  J) 

_ .c _  _ _ _ _ _ 

C  THIS  FUNCTION  CALCULATES  THE  BINOMIAL  COEFFICIENT  (  I  ) 

_ C _ _  _  _  (  J  ) 

C 

DIMENSION  BINOMI 15.15  ) 

. .  .  . . - . — . 

_ COEF  ■  1_. _ 

GO  TO  51 

_ 1  COEF  «  BINOMII.J) _ 

IFICOEF.NE.O.  )  RETURN 

i  F  ( J- 1  )  _2  * *3 _ 

2  COEF  ■  1 

_ GO  TO  31 _ 

A  CONTINUE  .  . 

_ I.rJ+ 1 _ 

COEF  ■  1. 

DO  10  L  =  K.I 

. . COE  >"  »"'co  FF*L . 

10  CONTINUE 
’DO  2C*N“Vl  ,  j 
COEF  ■  COEF/  N 

- ?-0  foN'riNur - - — - - 

BINOMI I.JJ  *  COEF 

- oPTOTTn - - - - -  -  - . . 

ENTRY  I N I T 

- - ira^rirriTJTs - - - 

DO  50  L*  1  » 15 

- griromTTrror-cr; — - - - - - 

30  CONTINUF 

- STTORTTUCfF - - -  '  '  - - 

RFTURN 

- FWT - - -  -  - . — 


h 

V. 


I 

f 


4  S; 

1 ■■ 

i 
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SUB ROU f I NE  f  NV£ RYfB  ,  K  ,  K  2  .  XMTN  , P L ACT 

C 

c 

this  subroutine  sets  up  a  unit  matrix  adjacent  to  b 

ELEMENTARY  ROW  OPERATIONS  ARE  THEN  PERFORMED  ON  THE 

(K.K) 

NEW 

K  X  2K  MATRIX 

c 

c 

TO  RFDUCE  B(K.K)  TO  A  UNIT  MATRIX.  THIS  WILL  PLACE 

the  matrix  bik.x)  in  the  right  half  of  bik.2K) 

THE 

INVERSE  OF 

c 

c 

c 

c 

ON  exit,  the  inverse  of  b  replaces  b 

6  IS  AN  ARRAY  OF  2 *K**2  LOCATIONS  CONTAINING  THE  MATRIX 

K  IS  THE  DIMENSION  OF  THE  SQUARE  MATRIX  B 

K2  IS  2*K 

c 

c 

XMIN  IS  THE  SMALLEST  ALLOWABLE  MAGNITUDE  OF  THE  PIVOT 
FLAG  WILL  BE  RETURNED  AS  0.  IF  THE  INVERSION  WENT  OFF  OK 

c 

c 

flag  will  be  returned  as  10.  if  a  pivot  element  was  too 
flag  should  be  tested  after  each  call  to  this  routine 

SMALL 

c 

DIMENSION  BIK.K2) 

c 

FLAG  =  0. 

c 

c 

SET  UP  UNIT  MATRIX 

c 

DO  1  1=1. K 

DO  1  J=1  .  K 

BII.K+J)  =  0. 

IF(I.EQ.J)  BII.K+JI  =  1. 

1  CONTINUE 

C 

c 

FIND  LEADING  ELEMENT  WITH  GREATEST  MAGNITUDE 

c 

DO  6  J=  1  »K 

M  =  J 

N  =  J  +  l 

DO  2  L=N.K 

IF  (ABS(BIM»Jn.LT.ABS(D(L.JIU  M  =  l 

C 

_C_ 

2  CONTINUE 

IF  (ABSIOIM.J) I.LT.XMIN)  GO  TO  10 

INTERCHANGE  JTH  AND  MTH  ROWS 

c 

_ D.O.J _ L- J  .  K2 _  .  _  .  . 

D  =  B(J.L) 

B(  J,L  )  =  8  <M  » L  ) 

DIM.L  )  =  0 

3  CONTINUE 

C 

C 

ZERO  OUT  PIVOTAL  JTH  COLUMN.  SKIPPING  PIVOTAL  JTH  ELEMENT 

C 

C 

DIVIDE  JTH  ROW  BY  PIVOT 

C 

DO  4  M=N»K2 

B(J.M)  =  B ( J  »M )  /  B(J.J) 

U  CONTINUE 

C 

DO  6  M=1  ,K 

c 

c 

M  DETERMINES  ROW  BEING  MODIFIED.  ONE  WHOLE  ROW  AT  A 

TIME 

IF  (  M.EQ.J  )  GO  TO  6 

DO  5  L=N.K2 

c 

C 

L  DETERMINES  ELEMENT  IN  THF  MTH  ROW 

c 

B(M.L)  =  B(M.L)  -  B(M.J)  *  BIJ.LI 

106 
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■>  CONTINUE 

_  6  CQNTJNUF 

C 

C _ INVERSE  OF  B  IS  NOW  IN  RIGHT  HALF  OF  BIK.K2) 

C  NOW  MOVE  B  INVERSE  TO  WHERE  B  WAS 
DO  .7  I  *  1  •  < 

DO  7  J=1  ,< 

B(  I  ,JI  ~  fill .J+K) 

7  CONTINUE 

_ RETURN 

10  FLAG  =  17. 

RETURN 
END  . 


107 


i 


FUNCTION  '■’INUS1  i  N  ) 

C 

"c"  THIS  FUNCTION  RF TURNS  <-ll»*N  ,  0.LE,N,LE.1R 

C _ 

DIMENSION  Ml  20) 

DATA  I (Mil). 1  =  1,20)=  1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 

r  i,-i  > 

MINUS  1  =  VKN  +  1  1 
RETURN 

_ gND . . 


I 

i 

l 


l 
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SUBROUTINE  NORMAL ( A .PROD .M.N I 

C~Ws"  SlTaROUT  Ine'  MUL  T I  PL  I ES  A*  (  A -TRANSPOSE )  *  PUTS  PRODUCT ‘INTO  PROD' 

C  A  IS  (MXN).  PROD  IS  THEREFORE  <M  X_  M) _  _ _ 

C  . 

DIMENSION  AIM.NI  .  PROU(M.M) 

c . .  . .  '  ‘  . .  . . 

DO  _1  I « 1 «  M 

do  l 

_ _ PR0DU..LI-0.  _ _ _ 

DO  1  J-l.N 

1  PROD( I  «L  )  _■  PROO(I.L)  +  A ( I »  J ) *A ( L  « J ) 

return" .  "  ‘  . . 

END 


*Su¥ROuTfNF”MXMULT_(AfB.'C.M.N.<  ) 

C 

C  THrs*SUBRbUT'lNE  MULTIPLIES  MATRIX  A  BY  MATRIX  B  AND  STORES  THE 
C  PRODUCT  IN  C.  (C  CANNOT  BE  THE  SAME  AS  A  OR  B. I 
C 

C  A  I S  (MX  N) 

r"B~Ts"(N~xXF 

C__  C_J_S._1M.X_K  ) 

c~ 

_ DIMENSION  A ( M «N )  .  B(N.KI  t  C(M,K)  _ 

C 

DO  1  1*1  .M 

WI  "l-TFk"  '  '  . . 

C( I»L)  = 

“"do  r  “=t.n 

C(  I  .L  )  *  C  (  I  »L  )  ♦  A(  1  ,J)  «  B  (  J  *  L.  >  _ _ _ _ __ 

1  CONTINUE 
RETURN 

----  . 


no 


t 


I 

p 

i 
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SUBROUTINE  MX TMUL ( A , B  .C *M.N  »K) 


DIMENSION  A(M.N).  B(M,K)  .  C(N.K) 

C  THIS  SUBROUTINE  MULTIPLIES  ( A-TRANSPOSE )  •  B»  PUTS  PRODUCT  IN  C 
C  A  IS  <M  X  N) 

_C_  A-TRANSPOSE  IS  _( N  X  M) 

‘c """if  7s  (m‘x  k7 

C  C  IS  (N  X  _K ) 

DO  1  I  *1  »N 


C(I,L)  *  0. 

DO  1  J»1.M 

TaT,rr-"cu VlT'+""ai j . n  *  bTj.u 

RETURN 


"SUBROUTINE  sYSxiTMBT ,11 . .“x.Y  .CONST) 

COMMON  /SIZES/  NEQ.  M.  MPR .  NX,  NUNK  .  NUNKP1 

“COMMON  "/TnI  T 7  INI  TA  V I  N  I  T  B  »  I  N I  TC  (  2  )  '»  I  N I  TO  ( 2_)  »  Tn"i"TK  1(2')  '."TfiTfKTfV 

_ 1 _ I_N  I  TK  3J  2J _ »  _I  N I T  K  4  .  IN1TK5.  INITK6I2)  .  INITK7, _ 

2  1nTtK8(2)""  INITK.9 InTtkIo,  TNfTKll  (2)  »  INITOMI2) 

COMMON  /l/  GI148.S71  .  L>N<56.S7>  .  E(2)  .  NUI2.1  . 

...  ------  -qy^-cTn*  VF‘ 

TYPE. REAL  NU 

DATA  €  IN  i  TA  =  I0!  5  !  I N  I  T B  =  1 6  )  ,  ( I N I TC=22 ,41 >  ,  (  INI  TD  =  29 , 48  )  . 

_ 1 _ (  1NITK1  -  1  .  3fi  )  .  (  I N  I  T  K  2  =  2.)  »  UN  ITK  3  =3  »39Jj>  (INITK4=4)  . 

2  ( INI TKS  =  5 )  .  ( IN  I T<6  =  6 . 40  )  ,  ( IN  I T  K  7  =  7 )  ,  ( I N I TK8  =  8 . 4 1  )  , 

_ 3  _ (  I N I TK9=9 )  ,  (INITK10=10>  _,  {  IN  I  T  K  1.1  =  1 1 ,42  )  .(  IN  I  TOM  =  37. 56  ) 

V 

C _ THE  VARIOUS  INITS  DEFINE  THE  POSITION  OF  THE  UNKNOWNS  IN  THE  G  MATRIX 

r 

C  INITA+1  WILL  BE  THE  COLUMN  OF_UNKNOWN_A1  .  REGION  1 _ 

C  INIT8+1  WILL  BE  THE  COLUMN  OF  UNKNOWN  Bl.‘ REGION  I 
C  IN  I  TC  (  J )  -*■  1  WILL  BE  THE  COLUMN  OF  Cl,  REGION  J.  J=I»II 

T  fNfTOrJ  r+T  W  ILL  'BE  ThE~  COLUMN  OF  01,  REGION  J.  J»T.Tl . . 

C 

c  TnoteThat  aV.b’Vci  uo  not  exist,  so  that  ai  over  lie  s  “TnT  t  RTT»  . 

c  B1  OVERLIES  A  7 ,  AND  Cl  OVERLIES  B7)  _ 

C  ‘  . . 

C  INITKJ  WILL  BE  column  OF  KJ  ,  J»2,4,5,7,9,10 

c  ‘TnTTkjTTT  wTll  Se  column' of  k j ,  region  l  ,  j  =  i,3,6,8,Ti''VX=T,Tr 

C  INITOM(J)  WILL  6E  COLUMN  OF  OMEGA  ZERO  IN  REGION  J  ,  J= I , I  I 

C  . .  .  .  . . 

C  A  ,  U  ♦  K2 .  K4,  Kb,  K7.  _K9,  JKlU,  DO  NOT  EX  1ST  IN  REGION  II _ 

c  . . 

G  (  I  ROW , IN  I TKH 1 1 )  )  =  G( IROW, INITK1  (  I  I  )  )+  CONST  *  6,*X*Y 

. G  (  I  ROW  ,Tn  I  T kTuT)  r  = '  G( IROW  ,  IN  l  TKS  (  I  I  )  )  +  CONST *  ?'• 

GO  TO  (5,10),  II 

p~fR0w7lNTfK2T . =  G<  IROW.  INI TK2  )  +  CONST*  2,*X . 

_ G  (  IROW.  INITK5) _ =  G(  IROW. INI TKS )  +  CONST*  fc.«Y _ 

Mil  =  MPR  . . 

_ GO  TO  .1.5 _ ....  _  _ 

TO  Mlf  V”m" 

_ IS  CO  AS  Ml  8  ]  ,  Mil 

M3  =  Ml"+  1 

_ DO  45  N2  =  1  .M3  _ _ _  _  _  _  _____ _ 

N1  =  N2  -  1 
I  F_(  Ml  .  NEjJN  1  )  GO  TO  20 

. ~"xd'  =  T." 

xc  =  x 

_ XA  =  XC _ _  _  _  __ _ 

GO  TO  . .  '"  . 

20  XD  =_X  *  *  (  2  *  (_M  1  -N 1  )  I 
XC~  >  *  XD 
xe  =  xd 
~XA"  =  "x“c 

2 S  IF(Nl.NE.O)  GO  TO  30 

YD  =1.  -  -  .  - . - .  -  - . 

YC  =0. 

y"b"  =y'c  .  . . . . 

YA  =YC 

-_0  T-0— . - 

30  IF(Nl.NE.l)  GO  TO  33 
YA  =  1. 

GO  TO  34 

33*  Y A "*  y5iT?*rVi-2) 

_ 34  YB  =  Y  *  YA 


112 


t 


» 


c 


YC  =  YR 
YD  «_  Y*  YC  _ 

Ys  G  (  IROwVlNITD(!I)  +M1  1  =  Gt  IROW,  INI  TD<  I  Il+Ml )+  CONST 

1  *  MINUSI (N1  )*2.*(N1+1  )*I2.*N1+1 )*  YD  *  XD 

2  *  COEFI2*Ml+r  »2*N1+1) 

IFIM1.EO.  1  )  GO  TO  44 


G(  IROW,  INITCI  I!  )+M1  )  =  GUROW*  IN  ITC  (  III +IY]  )  +  CONST 


»  MINUS1 (Ml )*2.*N1*< 2.  »N1  +  1 • ) 
*  COEF ( 2*M1+1  ,2*N]+1) 


GO  TO  <4  0,44  )•  I  I 

40  G(  IROW.  INI  TB+M1  )  =  G  (  I  ROW  ,  1  NI  TB+Ml) 

“T  *” MI NUS1  ("nTY*  2 •  *N  1  *  <  2 . *N 1 + 1  .  > 

2  *  C0EF(2*K1.2*N1  ) 

'GIT  ROwVi  ffi  T  A+M1)  V  GIIROW, INITA+M1  ) 
1  *  MINUS1 (Nl  )*2.*N1*(2.*N1-1 . ) 


2  *  COEFI  2*v,l  ,2*N1  ) 

44  CONTINUE 

"4  s  c  o  n  t"  i  nU  e 

RETURN 

end" . 


YC  *  XC 


CONST 
Y b  *  xb' 


CONST 
YA  *  XA 


subroutine  siGYriRow.fi  .xVyYcONST) 

_ TYPE  real  N1 _  _  .  _ _ _ 

""COMMON  /"sizes"/  "nEQ»"m»  MPR.  NX,'  NUNK  ’»  NUNKPl"" 

COMMON  / 1  /  0(146.57)  ._GN(  56.57)  .  E  ( 2  )  .  N'J  1  2  )  , _ 

1  AIN,  PIN.  CIN,  VF 

COMMON  /INI  1/  INITA.  INITB.INITC(2  )  » I N1 TD ( 2  )  ,  INIJKK2).  INITK2, 
”T  IN  ITK3  ( 2  )  .  I  N  I  T  X  4  .  IN1TK5."  INI  T  <6  (2  )  ,  TnTtkT,"" 

_ 2..  I  N  I  T<8  (  2  )  .  INITK9  ,  INITK10.  I  N  I  TK.1 1.(  2)...  INIJOMl  2  ) _ 

GUR5w,TNif<8(;i)  )  =  G<  TROW,  INI  TK8  (II)  !  +  CONST  *  6.*  X  *  Y 

G ( IROW. INI TK 11(11))  =  G (  I  ROW .INITKll(II)  ) +  CONST  ♦  2. _ 

GO  TO  (5  ,ioV,  fi 

1 5  _G(_IR0W,INITK9)^  =  G( I  ROW  ,  IN  I TK9 )  +  CONST  *  2.*  Y 

"gTIROwTI  N  iTkTu  )  =  CT I  ROW  »I  Nl  TK  10 )  +  CONST  *  6."*"x  . 

MU  =  MPR 

go"toTs 

10  Mil  =  M _ _  _  _ _ 

15  DO  45  Ml  =  1,  Mil 

__  ^3  =  Ml _ +  1 

00*45  N2  V  ]  ,M3 
N1  =  N2  -  1 

FfTmiTneVn  1 +  lT'GO  TO  20 

XD  =  1. _ _ _ _ _ _  _ _ 

XC  =  X  -  . - . 

XB_  =  XD 

xT“’xc 

GO  TO  30 

20"IFTmT. NEVnV) ' GOTO  25 
XD  =  0. 

XC  =  0.  "  ~ 

XB  =  0._ 

xTV'o.  . 

GO  TO  30 

Ts"  x'D"^"x  *“*  (  Ym'mT-n  I- 1 V ) 

XC  =  X»  XD _ _  _  _ _ 

XB  «  XD 
XA  =  XC 


30  IF(Nl.NE.O  GO  TO  3b 
YA  =  1. 


GO 

TO  40 

35 

YA 

=  Y  ** ( 2  *N  1  ) 

40 

YR 

=  Y  *  YA 

YC 

=  YB 

YD  =  Y  *  YB 

G< IROW.INITDJ I  I )+Ml  )  =  G ( ) ROW, IN  I T0( II ) +M1 ) +  CONST 

1  . TmTnUS'I  ( N 1)  *2  .  *  (  K 1  -N 1 )  «  ( 2  .  *  I M 1 -N 1  )  - 1.  )  *  YD'  *  XD  ” 

2  »  COES  ( 2*M I + 1  . 2*N1  I  )  _ 

I F  (  Ml .  El5*  1 )  GO  TO  44 

G(TroT,TnT  TC( f I) +M 1)  =  G {  I  ROW ,INITC(TI)+M1)+  CONST 

__  1  « ( 2»  *  I  M 1-N 1 )*I. )»?.»( Ml-Nl )»M  I NUS1 (N1 )  *  XC  *  YC 

.  coY(r  (  2 *M  1  f  ]  ,  2 *N  1  +  1  ) 

GO  TO  (42  ,44  ).  I  I 

42"G( IROW. INlTb+Ml )  =■  G( IHOa.INI Trf+Ml )  +  CONST 

1  *  MINUS1  C  Ml ) *2*(M1  — N1  ) * ( 2 . * ( M 1 — N 1 )-l . )  *YB  *  XB 

.  ---  __  coeF  (  2#M1  ,2*N1  ) 

G(  IROW,  IN  ITA+M1  )  =  C,(  IR0W.INITA+M1  )  +  CONST 

1  *  MINUSKN]  >*2.*(M1-N1)»(?.»(M1-\1  )  +  l.  )  *  XA  *  YA 

2  _ *  COEF(  2»M1 ,2*N1 ) 

44  CONTINUE 

45  CONTINUE 
. RETURN 

END 
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SURROUflNr"  I AUX Y (  I  ROW •  I  I  ,X ,Y  .CONST) 

COMMON  /\/  0(148,  57)  ,  (iNI  56.57)  ,  E  ( 2  )  ,  NU  ( 2  >  , 

T  '  A  1  N  *  B  1  N ,  CIN,  VF 
COMMON  /SllFl/  NFQ.  M,  MPK ,  NX,  NiJNK.  ,  NUNKP1 

COMMON  /  I N I T/  1NITA, iNlTH.INI  TC(P  )  ,  INI  TU(2)  »  Tn  I  TK1  <  2  )  ,  '  I  N  I  TIC2  . 

1  IN  I  TK? (2  >  ,  INITK4  ,  INITK5,  INIT<6(2>  ,  INITIO, 

2  IN  I  T<8 ( 2 )  ,  INITIO  ,  INITK10,  INITK1K2).  INITOM(2)~ 
TYPF  RFAL  NU 

C 

_ _ G ( IROW , IN  I T  K 1 ( 11))  =  G( I  ROW, INI TK1  (  II  )  )  - 3 .«  Y  *  «2_ «  CONST _ 

G( IROW. INITK3I I  I » )  =  G ( I  ROW • I N I TK 3 (  I  I))  -CONST 

G  (I  ROW  »  I N  I  TK8 ( I  I  I  I  =  G(  1  ROW,  INI  TICS  (III)  -  CONST  *  '4.*X**2 _ 

"GO “TO  7 5 , fJT.  ff" 

5  _G(  IROW,  INITX.2  >  =  G(  IROW,  IN  ITK2  )  -  2.*  CONST  *  Y 

'"gITrOwTTn I  T’k9)  =  G< IROW, INITK9)  -  2.*  CONST  *"x 

_ MU  =  MPp _  _ 

GO  TO  IS 
10  MU  =  M 
lT  DO  45  Ml”*  T,  M I  I 
M3  =  Ml  +  1 
DO  '45  N2~  =‘  1  ,  M3 
N1  *  N2  -  1 
IF(Ml.NE.Nl)  GO  TC  20 
XD  _  =  0. 

X  H  =  0. 

..... . . 

GO  TO  25 

20  XD  =  X**( 2*(M1~N1 1-1 1 
XC_f  X  _*  £D 

“Tb  *'75  ’  . . . . 

XA  *  XC 

r5"rr(Tr;NF"o"T"Go"fo  30 

_ YA  =  0. _ _ _ _ _ 

YB  =  1. 

_ GO  TO  35 _  _  .  _ _ _ _ _ 

30  YA  =  Y  **  (2*Ni-  if 

__  YB  =  Y  *_  YA  __  _ 

'*35"YC“="YB 

_ YD  =  Y»YC _ _ _ _ _ 

G( IROW, INI  TD< 1 1 I+Ml )  =  G<  IROW,  INITDU  I  l  +  Ml  I  -  CONST 
1  *  MINUSilNn  <M2.*Nl+2.  )  *  2.»(M1-N1>*  YD  *XD 

—j  *  _coefT2*mT + r ,  2  *N1+T  I  ‘  ' 

IFIMl.EQ.il  GO  TO  44 


G  (  I  ROW  ,  I N  1  TC  (  II  I  +M 1  )  *  G(  IROW,  INI  TC(  I  I  l+Ml  )  -  CONST  _ 

1  *  MINUSl(Nl)*  ( 2 • *N 1 +1 • )  * ( 2 . * ( Ml -N 1 1  +  1  .  )  *  YC  *  XC 

_ 2 _ »  COEF  ( _2*M  1  +  1,2 *N  1  + 1  I 

”  "go"to"  (‘40744  iVIT” 

40  G(  IROW.INITR+Ml)  *  G  (  I  ROW  ,  I N  I  TB+M1  )  -  CONST 

'I  MINUSlINlT*' (2.*N1  +  1  »*  2  .  *  (  Ml -N  if  *  yeT'*"X6‘ 

2  *  C0EF(2*M1 ,2*N1 I 

6( Ift6w,  IN  1  TA+M1 )  *  =  G< IROW, INlTA+MlT  -  CONST 

1  *  MINUS1 (N1 )*2.*N1*(2.« (Ml-Nl )+l. )*XA«YA 

-  ■*“c6efT2*MT;2‘*N'1  )  . '"r 

44  CONTINUE 

45  CONTTnUE  . 

RETURN _ _  _ _ _ 

END 
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SUBROUT  I Nfc  Ul  I  ROW* I  I .X.Y.CONST  ) 

COMMON  /l/  6(148,27)  ,  GN(Sf>,57)  ,  E ( 2 )  .  NU ( 2  )  . 

1  AIN.  BIN.  CIN,  VF 

COMMON  /SIZES/  NFO.  M.  MPR .  NX.  NUNK  ,  NUN K Pi 

Common  /init/  in i  t  a  « i  n  i  tb.T  n  i tc  1 2 1 » 1 ni to <2 > ,  1  Mi fX 172TT  in  1 TK2 » 

1  IN  I TK3  <  2  )  .  INITK4  .  INITK5.  INITK6I2)  .  INITK.7, 

2  ~ .  1 N 1 TK8  <  2 )  .  INITK9  ,  INITK10,  lNITKli(2).  INITOMIZ) 

TYPE  REAL  NU 

ONENU  =  1.  -  NUl  I  I) 

_ TWONU_=  2,  -  NUl  II)  _ _____ 

REVNU  =  0.  -  NUIIII 

r 

ECON= (  I  lT-1 7*(  ( 1  •  +  N U (  2  )  )  *E  ( 1  )  /  (  El  2  >  *  ( 1.+NIJI  1  )  )  )  -1  .  )+T.  )*  CONST 
C 

"GTTROW  .1  N  I  T  K 1  I  I  I  )  i  =  G(  I  ROW  •  I N  I  T  K  1  (II)  l+ECON*  (  3'.*ONENU*Y*X*'*2 
1  -  TW0NU*Y»*3  ) 

GiTROW  .  I‘nTTK3(TI  )  )  =  6(1  ROW. IN  I TK 3 (  I  I  )  ) +ECON»  I -TWONO»YT” . 

G ( I ROW , IN  1 TK6 ( I  I ) )  =  G ( I  ROW • IN  I TK6 ( II ) ) +FC ON* ( 2  •*ONENU*X ) 

■  G(  I  ROW-, INff<8(  II))  =  G(  I  ROW .  IN  T  T<8  (  I  I  )  >  +ECON*  (  3  .  *REVNU*Y#X**2 
1  -  ONENU  «  Y**3) 

'"6(  IROW.INITXIK  II  )  )=  G  (  IROW  .  I N  I  TK  1 1  C  II  )  )  +ECON*  (  2.*REVN(j*)C)” . . 

G(  I  ROVn.INI  TOM  (  I  I  )  )  =  GIIROW  .  I N  I  TOM  (  II)  )  +ECON*  (  ONENU*Y  ) 

- goto-  T^rrorm  ‘  - -  - - - 

3  G( IROW, INI  TK2  1  =  G<  1  ROW  ,  IN  I  TK2  )  +CCON*  (ONENU*X*  »2- TWONU*Y»«?  ) 

. GirROWVrNTTCA)  =  Gt  TROW,  INITK4  )+ECON*REVNU  '  ‘  - 

G< IROW, INITKS  )  =  G( I  ROW  » INI TK5 ) +ECON* (6 . *ONENU#X*Y ) 

•-  . . -GdROWVTNirKR')  =•  G  (  I  ROW  »  IN  I  T<9  )  +ECON*  (  ?  .  *REVNU*X*Y) - - - 

G(  I  ROW  .IN  I  TKl'l  )=  G  (  IROa’,  INI  TK10  )+ECON»(  -3.*0NENU*  Y»*? 

- 1 - >3.-*RFVNU«X**2)  -  - 

MU  =  MPR 

- 'GOTO-  1? .  .  — . . 

10  Mil  =  M 

TTDTT55 . Mr  *  1,  Mir  •  . * 

M3  =  Ml  +  1 

- TO--4r- . N?~=1 ,  -w? -  •  -  ^  ‘ 

Nl  =  N2  -  1 

f  PTmTI  nEVn  1 )  Go  To  ’  20 

_X02  =  0. 

XB2  =  C. 

XA2  I. 

GO  TO  22 

2 0  XD2  ‘  =''  X  " ##'(  ?#  (  u  1  -Nl  )  -  J  ) 

XC2  =  x  *  xr>2 
XB7  =  XT)2 
XA?  =  XC2 

- yTf1rr2-¥{;ol_r?1)+-1  ,  -  - 

XC1  =  X  *  XD1 

- X01  =~XD1  - - 

XA1  =  XC1 

IF ( Nl • NT . ' )  GO  TO  3N 
YD1  =  1. 

YC1  =  0. 

YDi  =  n. 

YA1  =  0. 

YA2  =  1. 

GO  TO  40 

OP  IKMl.NE.l  )  60  10  32 
Y01  =  Y**2 
Yf  1  =  Y 
YD1  -  v 
YA 1  =  1. 
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35 


YA2  = 
GO  TO 
Y  A1“ 
Y31  = 


Y**2 

AO 

Y**  (  2*N 1-2  ! 
Y*  YA1 


YC1  - 
YD1  = 
'  YA2  = 
AO  Y52  = 
‘  YC2  = 
Y02  = 


YB 1 

Y*  YC1 

VbT~ 

Y*YA2 

Yb2 

Y*yC2 


G(  I  ROW*  IMTDI  I  11+M1  )  =  o(  IROW.  IN  I  T  0  (  I  ll+Yl  )+  ECON 

1  *  MINUS  1  (  N1 )  *  <ONENU*< 2.+N1+1 . 1*2.* (Nl  +  1  . ) / <?.*(M1-N1  )+l. ) 

2  * 'Ybr«  XD1  ♦2.*RfeVNU  »  ( M 1  — m  1  )  *  XD2  *  YD2  ) 

3  *  CPEF(2*M1+1  ,2*N1+1) 


IF(M!.E0.1 )  GO  TO  AA 


G ( I  ROW  t IN  I TC ( II )+Ml)  =  G(IROW,!NITC(III+Mll+  ECON 

T .  minus l ( ni  )  *  <onenu*(2.*ni+i .>*  m l /  ( i*ii-ni+ i . )  "if" TCT'V'YiT 

2  +  REVNU  *  ( 2 •*  ( Ml-.N  1  )  +  1  •  )  *  YC2  »XC?  ) 

3  *  "COEF  ( 2*M  1  + 1  «2*N1  +1  ) 

GO  TO  (A?,AA)  .  I  I _  _  __ _ _ 

G < IROW *  I N t TB+M1 ) 


A 2  G ( I  ROW .  IN l TH+M1 )  =  G(  IROW. INITB+M1 )  +  ECON 

1  *  MINUS!  <N1  )*(ONENU*<  2.»N1+1  t*2.*Nl/(2.«  (Ml-Nl  )  + 1 .  )  * 

_  V“2'.  *REVN'  MM1-N1  I  »Y82  *  XB?  ) 

3  *  C0EF(2‘M1,2*N1) 

t-(TROW.'rNlYA+.Yi)  *  G< IROW.INITA+W1 !  +  ECON 
_ 1 _ «  MINUSHN1)  *  (O.NENU*  N1*<2.*N1-1>/  (Ml-Nl  +1.1 

2  '**  YA1*XSV  +'  'RtV~NU*<?.*<Ml-Nl>  +  l.) 

_3  «  COFF  (  ?*M1  «2*N  1  ) 

”a  A™  CONTINUE  . 

A 5  CONTINUE 


YH1«XB1 


*  YA2  *  XA2  T 


RETURN 

END 


/ 


117 


•>«tm  -  -- -ivf  ■  T  mr -  Wn  I  urn- 


SUOROl  IT  t  Nc  V  (  I  ROW,  !  I  .X,  v,  CONST  ) 

COMMON  /  1  /  G  (  148.57)  ,  GN(5A,S7)  .  F  (  2  >  ,  NIJI2)  , 

1  AIN.  BIN*  Cl  IX.  VF 

_ COMMON  /SIZFS/  N  FA-.  M.  VPR.  NX .  NUMC  .  NUNNPI 

C  GMf-’ON  "  /  I  N  1  T  /  IiJi7A.INITH.INI  TC<2  ),  INI  TD(  2  )  ,  IMITK1I2J.  IMTK2, 

1  INITK312)  .  INITA4  ,  INITK5,  i  N I  TK6 1  ?  )  .  INITY7. 

2  I N  I  T  X  8  ( ? )  .  INIfK.9  .  IMIKIfl.  IMITKIK2).  INITO.-K2I 
TYPF  REAL  NLI 

RFVNU  =  o.  -  NO (  I  I  I 
TWONIJ  =  2.  -  NU  (IP 
ONFNU  =  1.  -  NJ< I  I ) 

C 

“  tCON=  (  (  I  I  - 1  )  *  (  <  1  .+NLM2)  )  »F(  1  )/  1  F(  2  !  *  I  1  .+N'J<  1  ))  )-l  .  )  +  7  .  )*  CONST 
C 

. G<  f’ROW.  INITIO  (III)  =  C.  I  I  ROW-.  INITK3I  I  I  )  )  +tCON*Rr  VNIJ*X 

G ( I ROW. INI TK1  (III)  =  Of  I ROW. INITK1 ( I  I  ) J+FCON* <  3»RFVNU*X*Y**2 
1  -  ONFNU*  X**'l) 

G(IR0W.INIT<6(  II  )  )  =  G  ( I  ROW.  I N 1 1  X6  (  I  I  )  )+FCON*RFVNll*2.*Y 
"  G ( IRCW . IN  I TK 8 (III)  =  G( I  ROW . INITK8I I  I ) M-ECON* ( ONFNU*X*Y**2  *8. 

1  -  TWONIJ  *  X**3) 

‘SVfROWV  I N  I  TX  1 1  (  I  I)  )  =  G(  IROW.INITKUI  I  I  )  )  ♦  ECON»2.*ONENU* Y 
G  (  I  ROW  *  Ii\|  I  TOM  (II)'  =  G (  I  ROW  1 1 N I  TOM  (  I  I  )  )  +  ECON* ( -ONFNU*X > 

GO  TO  (5 , 1 U  j .  H 

_5  G  (  I  ROW  .INI  TK2  I  "  Gf  I  ROW,  IN  I  TO  )  +  ECON*  2.*RFVNU  *  X  *  Y 

GrfROwVi'NITV.S)=  f,(  I  ROW  »  I N I  T  X  A  )  +  t  CON*3*  (  -ONFNU*X*'*?+PE  VNU*  "Y**2") 
G  t I  ROW . IN  I TK7 )  =  P<  I  ROW, INITK7)  +  FCON*REVNU 

R"0w\  I N  I  TK 9  1  =  G  (  I  ROW » f  N I T  K  9 )  *  ECON*  <  ONFNU*  Y*«2-Ttf0N0*  X**21 
G ( I  ROW ♦ IN  I TK 1 0  )  =  G< I  ROW , IN  I TK IP >  +FCON*  ONENU*6.*X*Y 


-xrrr  s- 

mpr- 

GO  TO 

15 

10 

Mil  = 

M 

15 

DO  45 

Ml  =  1 

M '3" r  Ml'  ir~  I 

DO  AS  N2  =  1. 

Nl  '■  KiS"  —  l 

_IFIM1_.NE.N1!  GO  TO  20 

‘xaT"=  "61 

Xfll  =  0. 

'xcT~="o” 

XDl  =  0. 

XA2  =  X 
XB?  =  1. 

TfC?  ■  X" 

XD2  =  1. 

- WTFTfi - - 

20  IFIMl.NE.Nl+l >  GO  TO  25 
XF1  =  1. 

XDl  =  1. 

G6”TO”28 

25  XB1  =  X**<2* (Ml-Nl-1 I ) 

“TOT  -“xll 

28  XC1  =  X  *  XR1 
XA1  =  XC1 
XB2  _=__X*  X A 1 
X*A2  =”  X*  XB2 
XC2  •  X A 2 

. XD2*  XB?" 

20  IF(Y.NF.O)  GO  TO  35 
YA2  =0. 

YB2  =  1. 

YC2“*T.'  ' 

GO  TO  40 


t 


* 


A 


3s.  Y  A?  =  Y**(?*N1-1I 
v>3?  =  Y*YA2 
'  YC?  =  YB2 

40  YD?  =  Y»YR?  _  _ _  _ 

Y  A 1  =  YD? 

YH1.  =  Y*Y  A  i  _ _ 

YCl'  =  YB1 
YD1  =  Y*YC1 

G(  IROW.  INI  TD(  II  )  t  Ml  >  =  (j(  I  ROW  » I N I T  D !  I  i  )  +M 1 )  +  FCON 

1  »  MINUS1  (  NJ  >  *  <  ONr\U*2  .  *  (  Wi  1  -N  1  )  *  (  2.  *  ( ril.-Nl  ) -!.)/(  2.»Nl+3.  ) 

2  *  XD1  *  YD1  +  RE VNU*2 •  *(N1+I.)  »  YD2  *  XD2 ) 

3  _  *  COEF ( 2*M1+1  .2*Nl+]> 

IFlMl.EQ.l  )  GO  TO  44 

GUROW.INITCI  m+Ml)  =  G<  I  ROW,  IN  I  TCI  I  I  1+M1  1+  ECON 
"l  *~MINUS1(N1)  M0NFNU*(2.*(Y1-N1  )  +  l.  )*  ’(Ml-Nl  )7(N1+“iV) 

2  »  YC1  »  XC1  +  REVNU*  ( 2 •*N1  +  1_«  )  *  YC2  «XC2  )  _ 

3  *  C  0  E  F  ( 2  *  M 1  +  1  ,2*N1  +  1) 

GO  TO  (42 .44  )  *  II 

4?  G ( I  ROW  *  I  N  I  TB  +  Y1 )  =  G ( I  ROW  » I N I TB+M1  )  +  FCON 

1  *  MINUS1  INI  )*<ONENU*(Wl-\l  )  *  (  2  •  «  ( f'l )  -M  )-l  •  )  /(Nl+1.) 

..  V " x B 1  *  Yb  1  +  RFVNU  *  (?.*N1  +  1.)  *  Y 82  *  aB2) 

3  *  COEF(2*Yl  »2*N1 ) 

G  (  I  ROW  .INI  T>, +Y 1 )  =  “G  (  I  ROW  *  1  .si >  TA+M1  )  +  ECON  * 

1  *  MINUS1 (Nl)*(ONtNU*(2.*(Ml-M  )+l.)  *2.* (Ml-Nl ) / (  2.*N1  +  1 . > 

2  *~YA1"  *  XA1  +  RE VNU*2 •  *N  1  *  Y A2  *  XA2  ! 

;  *  COEF ( 2*M1 »2*N1 ) 

44  V  ON  TINUF 
4*i  CONTINI/F 
“RETURN 
END 
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